The underlying components that make up a computer can be simplified into three basic parts—the computing units, the memory units and the connections between them. In addition, each of these units has different properties that we can use to understand them. The computational unit has the property of how many computations it can do per second, the memory unit has the properties of how much data it can hold and how fast we can read and write to it, and finally the connections have the property of how fast they can move data from one place to the other.

Using these building blocks we can talk about a standard workstation at multiple levels of sophistication. For example, the standard workstation can be thought of as having a central processing unit (CPU) as the computational unit, connected to both the random access memory (RAM) and the hard drive as two separate memory units (each having different capacities and read/write speeds) and finally a bus that acts as the connections between all of these parts. However, we can also go into more detail and see that the CPU itself has several memory units in it, the L1, L2 and sometimes even the L3 and L4 cache, which have small capacities but very fast speeds (from several kilobytes to a dozen megabytes). These extra memory units are connected to the CPU with a special bus called the backside bus. Furthermore, new computer architectures generally come with new configurations (for example, Intel’s Nehalem CPU’s replaced the frontside bus with the Intel QuickPath Interconnect and restructured many connections). Finally, in both of these approximations of a workstation we have neglected the network connection, which is effectively a very slow connection to potentially many other computing and memory units!

To help untangle these various intricacies, let’s go over a brief description of these fundamental blocks.

Computing Units

The computing unit of a computer is the centerpiece of its usefulness—it provides the ability to transform any bits it receives into other bits or to change the state of the current process. CPU’s are the most commonly used computing unit, however graphics processing units (GPUs), which were typically used to speed up computer graphics but are not becoming more applicable for numerical applications, are gaining popularity due to their intrinsically parallel nature which allows many calculations to happen simultaneously. Regardless of its type, computing units take in a series of bits, for example bits representing numbers, and outputs another set of bits, for example representing the sum of those numbers. In addition to the basic arithmetic operations on integers and real numbers and bitwise operations on binary numbers, some computing units also provide very specialized operations, such as the “fused multiply add” operation which takes in three numbers, A,B,C, and returns the value A*B + C.

The main properties of interest in a computing unit are the number of operations it can do in one cycle and how many cycles it can do in one second. The first value is measured by its Instructions Per Cycle (IPC), while the latter value is measured by its clock speed. These two measures are always competing with each other when new computing units are being made. For example the Intel Core series has a very high IPC but a lower clock speed, while the Pentium 4 chip has the reverse. GPU’s, on the other hand, have a very high IPC and clock speed, however suffer from other problems which we will outline later.

Furthermore, while increasing clock speed almost immediately speeds up all programs running on that computational unit (because they are able to do more calculations per second), having higher IPC’s can also drastically effect computing by changing the level of vectorization that is possible. Vectorization is when a CPU is provided with multiple pieces of data at a time and is able to operate on all of them at once. This sort of CPU instruction is known as SIMD (Single Instruction, Multiple Data).

In general, computing units have been advancing quite slowly in the past years (see Figure 1-1). Clock speeds and IPC have both been stagnant because of the physical limitations of making transistors smaller and smaller. As a result, chip manufactures have been relying on other methods to gain more speed, including hyperthreading, more clever out-of-order execution and multi-core architectures.

Hyperthreading presents a virtual second CPU to the host operating system and clever hardware logic tries to interleave two threads of instructions into the execution units on a single CPU. When successful gains of up to 30% over a single thread can be achieved. Typically this works well when the units of work across both threads use different types of execution unit - e.g. one performs floating point operations and the other performs integer operations.

Out-of-order execution enables a compiler to spot that some parts of a linear program sequence do not depend on the results of a previous piece of work - therefore both pieces of work could potentially occur in any order or at the same time. As long as sequential results are presented at the right time, the program continues to execute correctly, even though pieces of work are computed out of their programmed order. This enables instructions to execute when others might be blocked (e.g. waiting for a memory access), allowing greater overall utilisation of the available resources.

Finally, and most important for the higher-level programmer, is the prevalence of multi-core architectures. These architectures include multiple CPU units within the same unit which increases the total capability without running into barriers in making each individual unit faster. This is why it is currently hard to find any machine with less than two cores—in this case the computer has two physical computing units which are connected to each other. While this increases the total number of operations that can be done per second, it introduces intricacies in fully utilizing both computing units at the same time.

Simply adding more cores to a CPU does not always speed up a program’s execution time. This is because of something known as Amdahl’s law. Simply stated, Amdahl’s law says that if a program designed to run on multiple cores has some routines that must run on one core, this will be the bottleneck for the final speedup by allocating more cores.

For example, if we had a survey we wanted 100 people to fill out, and that survey took 1 minute to complete, we could complete this task in 100 minutes if we had one person asking the questions (this person goes to participant 1, asks the questions, waits for the responses, then moves to participant 2). This method of having one person asking the questions and waiting for responses is similar to a serial process. In serial processes, we have operations being satisfied one at a time, each one waiting for the previous operation to complete.

However, we could perform the survey in parallel if we had two people asking the questions which would let us finish the process in only 50 minutes. This can be done because each individual person asking the questions does not need to know anything about the other people asking questions. As a result, the task can be easily split up without having any dependency between the question askers.

The speedup of adding more people asking the questions will give us speedups until we have 100 people asking questions. At this point, the process would take 1 minute and is simply limited by the time it takes the participant to answer questions. Adding more people asking questions cannot gain any further speedups because these extra people will have no tasks to perform—all the participants are already being asked questions! At this point, the only way to reduce the overall time to run the survey is to reduce the amount of time it takes for an individual survey, the serial portion of the problem, to complete. Similarly with CPU’s, we can add more cores that can perform various chunks of the computation as necessary until we reach a point where the bottleneck is a specific core finishing its task. In other words, the bottleneck in any parallel calculation is always the smaller serial tasks that are being spread out.

Furthermore, a major hurdle with utilizing multiple cores in Python is Python’s use of a Global Interpreter Lock (or GIL). The GIL makes sure that a Python process can only run one instruction at a time, regardless of the number of cores it is currently using. This means that even though some Python code has access to multiple cores at a time, only one core is running a Python instruction at any given time. Using the previous example of a survey, this would mean that even if we had 100 question askers, only one could ask a question and listen to a response at a time. This effectively removes any sort of benefit from having multiple question askers! While this may seem like quite a hurdle, especially if the current trend in computing is to have multiple computing units rather than having faster ones, this problem can be avoided by using other standard library tools like multiprocessing or technologies such as numexpr, cython or distributed models of computing.

Figure 1-1. Clock speed of CPU’s over time (Data from CPU DB)

Memory Units

Memory units in computers are used to store bits. This could be bits representing variables in your program, or bits representing the pixels of an image. In this way, the abstraction of a memory unit applies to the registers in your motherboard as well as your RAM and hard drive. The one major difference between all of these types of memory units is the speed at which they can read/write data. To make things more complicated, the read/write speed is heavily dependent on the way that data is being read.

For example, most memory units perform much better when they read one large chunk of data as opposed to many small chunks (this is referred to as sequential read versus random data). If data in these memory units are thought of like pages in a large book, this means that most memory units have better read/write speeds when going through the book page by page rather than constantly flipping from one random page to another. While this fact is generally true across all memory units, the amount that this affects each type is drastically different.

In addition to the read/write speeds, memory units also have latency which can be characterized as the time it takes the device to find the data that is being used. For a spinning hard drive, this latency can be high because the disk needs to physically spin up to speed and the read head must move to the right position. On the other hand, for RAM this can be quite small because everything is solid state. Below is a short description of the various memory units that are commonly found inside a standard workstation in order of read/write speeds,

• Spinning Hard Drive - Long term storage that persists even when the computer is shut down. Generally has slow read/write speeds because the disk must be physically spun and moved. Degraded performance with random access patterns but very large capacity (terabyte range)
• Solid State Hard Drive - Similar to a spinning hard drive however with faster read/write speeds but smaller capacity (gigabyte range)
• RAM - Used to store application code and data (such as any variables being used). It has fast read/write characteristics and performs well with random access patterns, however is generally limited in capacity (gigabyte range)
• L1/L2 Cache - Extremely fast read/write write speeds. Data going to the CPU must go through here. Very small capacity (kilobyte range)

Figure 1-2. Characteristic values for different types of memory units (values from Febuary, 2014)

A very immediate trend can be seen in the above list that read/write speeds and capacity are inversely proportional—as we try to increase speed capacity gets reduced. Because of this, many systems implement a tiered approach to memory: data starts in its full state in the hard drive, part of it moves to the RAM then a much smaller subset to the L1/L2 cache. This method of tiering enables programs to keep memory in different places depending on access speed requirements. When trying to optimize the memory patterns of a program, we are simply optimizing which data is placed where, how it is laid out (in order to increase the number of sequential reads) and how many times it is moved between the various locations. In addition, methods such as asynchronous I/O and preemptive caching provide ways to make sure that data is always where it needs to be without having to waste computing time—most of these processes can happen independently while other calculations are being performed!

Communications Layers

Finally, let’s look at how all of these fundamental blocks communicate with each other. There are many different modes of communication, but they are all variants on a thing called a bus.

The frontside bus, for example, is the connection between the ram and the L1/L2 cache. Here, we move data that is ready to be transformed by the processor into the staging ground to get ready for calculation and move finished calculations out. On the other hand, there are other buses such as the external bus which acts as the main route from hardware devices (such as hard drives and networking cards) to the CPU and system memory. This bus is generally slower than the frontside.

In fact, many of the benefits to the L1/L2 cache come from the fact that they are on a faster bus. By being able to queue up data necessary for computation in large chunks on a slow bus (from RAM to cache) then having it available at very fast speeds from the backside bus (from cache to CPU), the CPU is able to do more calculations without waiting such a long time.

Similarly, many of the drawbacks of using a GPU come from the bus it is connected on: since the GPU is generally a peripheral device, it communicates through the PCI bus which is much slower than the frontside bus. As a result, getting data in and out of the GPU can be quite a taxing operation. In this way, the advent of heterogeneous computing, or computing blocks that have both a CPU and a GPU on the frontside bus, aim at reducing the data transfer cost and making GPU computing more of an available option, even when a lot of data must be transferred.

In addition to the communications blocks within the computer, the network can be thought of as yet another communication block. This block, however, is much more pliable that the ones discussed previously; a network device can be connected to a memory device, such as a network attached storage (NAS), or another computing block, as in a computing node in a cluster. However, network communications are generally much slower than the other types of communications than were mentioned before. While the frontside bus can transfer dozens of gigabits per second, the network is limited to the order of several dozen megabits.

It is clear then that the main property of a bus is its speed: how much data can it move in a given amount of time. This property is given by combining several quantities: how much data can be moved in one transfer (bus width) and how many transfers it can do per second (bus frequency). It is important to note that the data moved in one transfer is always sequential: a chunk of data is read off of the memory and moved to a different place. Thus, the speed of a bus is broken up into these two quantities because individually they can affect different aspects of computation: a large bus width can help vectorized code (or any code that sequentially reads through memory) by being able to move all the relevant data in one transfer.

On the other hand, having a small bus width but a very high frequency of transfers can help code that must do many reads from random parts of memory. Interestingly, one of the ways that these properties are changed by computer designers is by the physical layout of the motherboard: when chips are placed closed to one another, the length of the physical wires joining them is smaller and can allow for faster transfer speeds. In addition, the number of wires themselves dictate the width of the bus (giving real physical meaning to the term!).

Figure 1-3. Connection speeds of various common interfaces (Made by Leadbuffalo, [CC BY-SA 3.0])

Understanding the basic components of a computer is not enough to fully understand the problems of performance programming. The interplay of all of these components and how they work together to solve a problem introduces the extra levels of complexity that are necessary. In this section, we will try to explore some toy problems and understand how the ideal solution would work and how Python approaches it.

A warning: this section may seem bleak—most of the remarks seem to say that Python is natively incapable of dealing with the problems of performance. This is untrue for two reasons. Firstly, in all of these “components of performant computing” we have neglected one very important component: the developer. What native Python may lack in performance it gets back right away with speed of development. Furthermore, throughout the book we will introduce modules and philosophies that can help mitigate many of the problems described below with relative ease. With both of these two aspects combined, we will keep the fast development mindset of Python while removing many of the performance constraints.

import math
def check_prime(number):
    sqrt_number = math.sqrt(number)
    number_float = float(number)
    for i in xrange(2, int(sqrt_number)+1):
        if (number_float / i).is_integer():
            return False
    return True

print "check_prime(10000000) = ", check_prime(10000000) # False
print "check_prime(10000019) = ", check_prime(10000019) # True

import math
def check_prime(number):
    sqrt_number = math.sqrt(number)
    number_float = float(number)
    numbers = range(2, int(sqrt_number)+1)
    for i in xrange(0, len(numbers), 5):
        result = (number_float / numbers[i:(i+5)]).is_integer() #not valid Python
        if any(result):
            return False
    return True

def search_fast(haystack, needle):
    for item in haystack:
        if item == needle:
            return True
    return False

def search_slow(haystack, needle):
    return_value = False
    for item in haystack:
        if item == needle:
            return_value = True
    return return_value

Python is highly expressive and easy to learn - new programmers quickly discover that they can do quite a lot in a short space of time. Many Python libraries wrap tools written in other languages to make it easy to call other systems, for example the scikit-learn machine learning system wraps LIBLINEAR and LIBSVM (both of which are written in C) and the numpy library includes BLAS and other C and Fortran libraries. As a result, python code that properly utilizes these modules can indeed be as fast as comparable C code.

Python is described as “Batteries Included” as many important and stable libraries are built-in:

Outside of the core language there are a huge variety of libraries:

A large selection of managed environments and shells are available to fit various deployment scenarios:

One of Python’s main strengths is that it enables fast prototyping on an idea. Due to the wide variety of supporting libraries it is easy to test if an idea is feasible, even though the first implementation might be rather flaky.

If you want to make your mathematical routines faster then look to numpy. If you want to experiment with machine learning then try scikit-learn. If you are cleaning and manipulating data then pandas is a good choice.

In general, it is sensible to raise the question “if our system runs faster - will we as a team run slower in the long-run?”. It is always possible to squeeze more performance out of a system if enough man-hours are invested, this might lead to brittle and poorly-understood optimizations that ultimately trip the team up.

One example might be the introduction of Cython (Cython), a compiler-based approach to annotating Python code with C-like types so the transformed code can be compiled using a C compiler. While the speed gains can be impressive (often achieving C-like speeds with relatively little effort), the cost of supporting this code will increase. In particular it might be harder to support this new module as team members will need a certain maturity in their programming ability to understand some of the trade-offs that have occurred when leaving the Python virtual machine that introduced the performance increase.

The first aim of profiling is to test a representative system to identify what’s slow (or using too much RAM, disk I/O or network I/O). Profiling typically adds an overhead (10* to 100* slowdowns can be typical) and you still want your code to be used as similarly to a real world situation as possible. Extract a test case and isolate the piece of the system that you need to test, preferably it’ll have been written to be in its own set of modules already.

The basic techniques that are introduced first in this chapter include the %timeit magic in IPython, time.time() and a timing decorator. You can use these techniques to understand the behavior of statements and functions.

We will cover cProfile (Using the cProfile module) to show you how to use this built-in tool to understand which functions in your code take the longest to run. This will give you a high level view of the problem so you can direct your attention to the critical functions.

Next we’ll look at line_profiler (line_profiler for line-by-line measurements) which will profile your chosen functions on a line-by-line basis. The result will include a count of the number of times each line is called and the percentage of time spent on each line. This is exactly the information you need to understand what’s running slowly and why.

Armed with the results of line_profiler you’ll have the information you need to move on to using a compiler (Chapter 7).

perf stat is later introduced in Example 6-8, there you’ll learn how to understand the number of instructions that are ultimately executed on a CPU and how efficiently the CPU’s caches are utilized. This allows for advanced-level tuning of matrix operations. You should read this when you’re done with this chapter.

After line_profiler we show you heapy (Inspecting objects on the heap with heapy) which can track all of the objects inside Python’s memory - this is great for hunting strange memory-leaks. If you’re working with long-running systems then dowser (Dowser for live graphing of instantiated variables) will interest you - you can introspect live objects in a long-running process via a web-browser interface.

To understand why your RAM usage is high we’ll show you memory_profiler (memory_profiler for diagnosing memory usage). It is particularly useful to chart RAM usage over time on a labelled chart to explain to colleagues why certain functions use more RAM than expected.

Finally we’ll give you an introduction to the Python bytecode inside CPython (The dis module to examine CPython bytecode) so you can understand what’s happening “under the hood”. In particular having an understanding of how Python’s stack based virtual machine operates will help you to understand why certain coding styles run more slowly than others.

Before the end we review how to integrate unit tests whilst profiling (Unit testing during optimization to maintain correctness) to preserve the correctness of your code whilst you make it run more efficiently.

We’ll finish with a discussion of profiling strategies (Strategies to profile your code successfully) so you can reliably profile your code and gather the correct data to test your hypotheses. Here you’ll learn about how dynamic CPU frequency scaling and features like TurboBoost can skew your profiling results and how they can be disabled.

To walk through all of these steps we need an easy-to-analyze function. The next section introduces the Julia Set - it is a CPU-bound function that’s a little hungry for RAM, it also exhibits non-linear behavior (so we can’t easily predict the outcomes) which means we need to profile it at run-time rather than analyze it off-line.

The Julia Set ^[4] is an interesting CPU-bound problem for us to begin with, it is a fractal sequence that generates a complex output image. It is named after Gaston Julia.

The code that follows is a little longer than a version you might write yourself, it has a CPU-bound component and a very explicit set of inputs. This configuration allows us to profile both the CPU usage and the RAM usage so we understand which parts of our code are consuming two of our scarce computing resources. This implementation is deliberately sub-optimal so we can identify memory-consuming operations and slow statements. Later in this chapter we fix a slow logic statement and a memory-consuming statement and then in Chapter 7 we significantly speed-up the overall execution time of this function.

We will analyse a block of code that produces both a false greyscale plot (Figure 2-1) and a pure greyscale variant of the Julia Set (Figure 2-3), at the complex point c=-0.62772-0.42193j. A Julia Set is produced by calculating each pixel in isolation, this is an “embarrassingly parallel problem” as no data is shared between points.

If we chose a different c then we’d get a different image. The location we have chosen has regions that are quick to calculate and others that are slow to calculate, this is useful for our analysis.

The problem is interesting because each pixel is calculated by applying a loop that could be applied an indeterminate number of times. On each iteration we test to see if this coordinate’s value escapes towards infinity or if it seems to be held by an attractor. Coordinates that cause few iterations are colored darkly, a high number of iterations are colored white. White regions are more complex to calculate and so take longer to generate.

Figure 2-1. Julia Set plot with a false greyscale to highlight detail.

We define a set of z coordinates that we’ll test. The function that we calculate squares the complex number z and adds c:

We iterate on this function whilst testing to see if the escape condition holds using abs. If the escape function is False then we break out of the loop and record the number of iterations we performed at this coordinate. If the escape function is never False then we stop after maxiter iterations. We will later turn this z’s result into a coloured pixel representing this complex location.

In pseudocode it might look like:

for z in coordinates:
    for iteration in range(maxiter):  # limited iterations per point
        if abs(z) < 2.0:  # has the escape condition been broken?
            z = z*z + c
        else:
            break
    # store the iteration count for each z and draw later

To explain this function, let’s try two coordinates.

Let us use the co-ordinate that we draw in the top-left corner of the plot at -1.8-1.8j. We must test abs(z) < 2 before we can try the update rule:

z = -1.8-1.8j
print abs(z)

2.54558441227

We can see that for the top-left co-ordinate the abs(z) test will be False on the 0th iteration, so we do not perform the update rule. The output value for this co-ordinate is 0.

Now let’s jump to the center of the plot at z=0+0j and try a few iterations:

c = -0.62772-0.42193j
z = 0+0j
for n in range(9):
    z = z*z + c
    print "{}: z={:33}, abs(z)={:0.2f}, c={}".format(n, z, abs(z), c)

0: z=              (-0.62772-0.42193j), abs(z)=0.76, c=(-0.62772-0.42193j)
1: z=    (-0.4117125265+0.1077777992j), abs(z)=0.43, c=(-0.62772-0.42193j)
2: z=(-0.469828849523-0.510676940018j), abs(z)=0.69, c=(-0.62772-0.42193j)
3: z=(-0.667771789222+0.057931518414j), abs(z)=0.67, c=(-0.62772-0.42193j)
4: z=(-0.185156898345-0.499300067407j), abs(z)=0.53, c=(-0.62772-0.42193j)
5: z=(-0.842737480308-0.237032296351j), abs(z)=0.88, c=(-0.62772-0.42193j)
6: z=(0.026302151203-0.0224179996428j), abs(z)=0.03, c=(-0.62772-0.42193j)
7: z= (-0.62753076355-0.423109283233j), abs(z)=0.76, c=(-0.62772-0.42193j)
8: z=(-0.412946606356+0.109098183144j), abs(z)=0.43, c=(-0.62772-0.42193j)

We can see that each update to z for these first iterations leaves it with a value where abs(z) < 2 is True. For this co-ordinate we can iterate 300 times and still the test will be True. We cannot tell how many iterations we must perform before the condition becomes False (and possibly this is an infinite sequence) so the maximum iteration (maxiter) break clause will stop us iterating potentially forever.

In Figure 2-2 we see the first 50 iterations of the above sequence. For 0+0j (solid line with circle markers) the sequence appears to repeat every 8th iteration but each sequence of 7 calculations has a minor deviation to the previous sequence - we can’t tell if this point will iterate forever within the boundary condition or for a long time or maybe for just a few more iterations. The dashed cutoff line shows the boundary at +2.

Figure 2-2. Two co-ordinate examples evolving for the Julia set.

For -0.82+0j (the dashed line with diamond markers) we can see that after the 9th update the absolute result has exceeded the +2 cutoff and we’d stop updating this value.

Below we break down the code that generates the Julia Set, we’ll analyse it in various ways throughout this chapter. As shown in Example 2-1 at the start of our module we import the time module for our first profiling approach and define some coordinate constants.

"""Julia set generator without optional PIL-based image drawing"""
import time

# area of complex space to investigate
x1, x2, y1, y2 = -1.8, 1.8, -1.8, 1.8
c_real, c_imag = -0.62772, -.42193

To generate the plot we create two lists of input data - the first is zs (complex z co-ordinates), the second is cs (a complex initial condition). Neither list varies and we could optimise cs to a single c value as a constant. The rationale for building two input lists is that we have some reasonable-looking data to profile when we profile RAM usage later in this chapter.

To build the zs and cs lists we need to know the coordinates for each z, in Example 2-2 we build this up using xcoord and ycoord and a specified x_step and y_step. The somewhat verbose nature of this setup is useful when porting the code to other tools (e.g. to numpy) and to other Python environments as it helps to have everything very clearly defined for debugging.

def calc_pure_python(desired_width, max_iterations):
    """Create a list of complex co-ordinates (zs) and complex
    parameters (cs), build Julia set and display"""
    x_step = (float(x2 - x1) / float(desired_width))
    y_step = (float(y1 - y2) / float(desired_width))
    x = []
    y = []
    ycoord = y2
    while ycoord > y1:
        y.append(ycoord)
        ycoord += y_step
    xcoord = x1
    while xcoord < x2:
        x.append(xcoord)
        xcoord += x_step
    # build a list of co-ordinates and the initial condition for each cell.
    # Note that our initial condition is a constant and could easily be removed,
    # we use it to simulate a real-world scenario with several inputs to our function
    zs = []
    cs = []
    for ycoord in y:
        for xcoord in x:
            zs.append(complex(xcoord, ycoord))
            cs.append(complex(c_real, c_imag))

    print "Length of x:", len(x)
    print "Total elements:", len(zs)
    start_time = time.time()
    output = calculate_z_serial_purepython(max_iterations, zs, cs)
    end_time = time.time()
    secs = end_time - start_time
    print calculate_z_serial_purepython.func_name + " took", secs, "seconds"

    # this sum is expected for 1000^2 grid with 300 iterations
    # it catches minor errors we might introduce when we're
    # working on a fixed set of inputs
    assert sum(output) == 33219980

Having built the zs and cs lists we output some information about the size of the lists and calculate the output list via calculate_z_serial_purepython. Finally we sum the contents of output and assert that it matches the expected output value. Ian uses it here to confirm that no errors creep into the book.

As the code is deterministic we can verify that the function works as we expect by summing all the calculated values. This is useful as a sanity check because when we make changes to numerical code it is very sensible to check that we haven’t broken the algorithm. Ideally we would use unit tests and we’d test more than one configuration of the problem.

Next in Example 2-3 we define the calculate_z_serial_purepython function, it expands on the algorithm we discussed above. Notably we also define an output list at the start which has the same length as the input zs and cs lists. You may also wonder why we’re using range rather than xrange - this is so in memory_profiler for diagnosing memory usage we can show how wasteful range can be!

def calculate_z_serial_purepython(maxiter, zs, cs):
    """Calculate output list using Julia update rule"""
    output = [0] * len(zs)
    for i in range(len(zs)):
        n = 0
        z = zs[i]
        c = cs[i]
        while abs(z) < 2 and n < maxiter:
            z = z * z + c
            n += 1
        output[i] = n
    return output

Now we call the calculation routine in Example 2-4, by wrapping it in a __main__ check we can safely import the module without starting the calculations for some of the profiling methods. Note that here we’re not showing the method used to plot the output.

if __name__ == "__main__":
    # Calculate the Julia set using a pure Python solution with
    # reasonable defaults for a laptop
    calc_pure_python(desired_width=1000, max_iterations=300)

Once we run the code we see some output about the complexity of the problem:

# running the above produces:
Length of x: 1000
Total elements: 1000000
calculate_z_serial_purepython took 12.3479790688 seconds

In the false-color plot (Figure 2-1) the high-contrast color changes give us an idea of where the cost function is slow-changing or fast-changing. Here in Figure 2-3 we have a linear color map, black is quick to calculate and white is expensive to calculate.

Figure 2-3. Julia plot example using a pure greyscale.

By showing two representations of the same data we can see that lots of detail is lost in the linear mapping. Sometimes it can be useful to have various representations in mind when investigating the cost of a function.

At the end of the code example above we see the output that is generated by several print statements. On Ian’s laptop this code takes approximately 12 seconds to run using CPython 2.7. It is useful to note that there is always some variation in execution time. You must observe the normal variation when you’re timing your code or you might incorrectly attribute an improvement in your code simply to a random variation in execution time.

Your computer will be running other tasks such as accessing the network, disk or RAM and these factors can cause variations in the execution time of your program.

Ian’s laptop is a Dell E6420 with an Intel Core I7-2720QM CPU (2.20GHz, 6MB cache, Quad Core) and 8GB of RAM running Ubuntu 13.10.

In calc_pure_python we can see several print statements. This is the simplest way to measure the execution time of a piece of code inside a function. It is a very basic approach, despite being quick and dirty it can be very useful when you’re first looking at a piece of code.

Using print statements is commonplace when debugging and profiling code. It quickly becomes unmanageable but is useful for short investigations. Try to tidy them up when you’re done with them or they will clutter your stdout.

A slightly cleaner approach is to use a decorator - here we add 1 line of code above the function that we care about. Our decorator can be very simple and just replicate the effect of the print statements. Later we can make it more advanced.

In Example 2-5 we define a new function timefn which takes a function as an argument, the inner function measure_time takes *args (a variable number of positional arguments) and **kwargs (a variable number of key/value arguments) and passes them through to fn for execution. Around the execution of fn we capture time.time() and then print the result along with fn.func_name. The overhead of using this decorator is small but if you’re calling fn millions of times the overhead might become noticeable.

We use @wraps(fn) to expose the function name and docstring to the caller of the decorated function (otherwise we would see the function name and docstring for the decorator, not the function it decorates).

from functools import wraps

def timefn(fn):
    @wraps(fn)
    def measure_time(*args, **kwargs):
        t1 = time.time()
        result = fn(*args, **kwargs)
        t2 = time.time()
        print ("@timefn:" + fn.func_name + " took " + str(t2 - t1) + " seconds")
        return result
    return measure_time


@timefn
def calculate_z_serial_purepython(maxiter, zs, cs):
    ...

When we run the above version (we keep the print statements from before) we can see that the execution time in the decorated version is very slightly quicker than the call from calc_pure_python, this is due to the overhead of calling a function (the difference is very tiny).

Length of x: 1000
Total elements: 1000000
@timefn:calculate_z_serial_purepython took 12.2218790054 seconds
calculate_z_serial_purepython took 12.2219250043 seconds

We can use the timeit module as another way to get a coarse measurement of the execution speed of our function. More typically you would use this when timing different types of simple expressions as you experiment with ways to solve a problem. Here we can use it to time the execution time for our CPU-bound function.

From the command line you can run timeit using:

$ python -m timeit -n 5 -r 5 -s "import julia1" "julia1.calc_pure_python( desired_width=1000, max_iterations=300)"

and note that you have to import the module as a setup step using -s as calc_pure_python is inside that module. timeit has some sensible defaults for short sections of code, for longer-running functions it can be sensible to specify the number of loops which are averaged (-n 5) and a number of repetitions of the loops (-r 5). The best result of all the repetitions is given as the answer, each result is an average of the loops in that test.

By default if I run timeit without specifying -n and -r then it runs 10 loops with 5 repetitions and this takes 6 minutes, so overriding the defaults can make sense if you want to get your results a little faster.

We’re only interested in the best case results as the average and worst case are probably a result of interference by other processes. Choosing the best of five repetitions of five averaged results should give us a fairly stable result. Try running the benchmark several times to check if you get varying results - you may need more repetitions to settle on a stable fastest-result time. There is no “correct” configuration - if you see a wide variation in your timing results then do more repetitions until your final result is stable.

5 loops, best of 5: 13.1 sec per loop

The overall cost of calling calc_pure_python is 13.1 seconds (as the best case) whilst single calls to calculate_z_serial_purepython take 12.2 seconds as measured by the @timefn decorator. The difference is mainly the time taken to create the lists of zs and cs.

Inside IPython we can use the magic %timeit in the same way. If you are developing your code interactively in IPython and the functions are in the local namespace (probably because you’re using %run) then you can use:

%timeit calc_pure_python(desired_width=1000, max_iterations=300)

It is worth considering the variation in load that you get on a normal computer. Many background tasks are running (e.g. Dropbox, backups) that could impact the CPU and disk resources at random. Scripts in web pages can also cause unpredictable resource usage. Figure 2-4 shows the single CPU being used at 100% for some of the timing steps above, the other cores on this machine are each lightly working on other tasks.

Figure 2-4. System Monitor on Ubuntu showing variation in background CPU usage whilst we time our function.

Occasionally the System Monitor shows spikes of activity on this machine, it is sensible to watch your System Monitor to check that nothing else is interfering with your critical resources (CPU, Disk, Network).

We can step outside of Python for a moment to use a standard system utility on Unix-like systems. The following will record various views on the execution time of your program and it won’t care about the internal structure of your code. Note that we specifically use /usr/bin/time rather than time so we get the system’s time and not the simpler (and less useful) version built into our shell.

If you try time --verbose and you get an error, you’re probably looking at the shell’s built-in time command and not the system command.

$ /usr/bin/time -p python julia1_nopil.py
Length of x: 1000
Total elements: 1000000
calculate_z_serial_purepython took 12.7298331261 seconds
real 13.46
user 13.40
sys 0.04

Using the -p portability flag we get three results:

By adding user and sys you get a sense of how much time was spent in the CPU. The difference between this and real might tell you about the amount of time waiting for I/O, it might also suggest that your system is busy running other tasks that are distorting your measurements.

time is useful because it isn’t specific to Python. It includes the time taken starting the python executable which might be significant if you start lots of fresh processes (rather than having a long-running single process). If you often have short-running scripts where the start-up time is a significant part of the overall runtime then time can be a more useful measure.

We can add the --verbose flag to get even more output:

$ /usr/bin/time --verbose python julia1_nopil.py
Length of x: 1000
Total elements: 1000000
calculate_z_serial_purepython took 12.3145110607 seconds
        Command being timed: "python julia1_nopil.py"
        User time (seconds): 13.46
        System time (seconds): 0.05
        Percent of CPU this job got: 99%
        Elapsed (wall clock) time (h:mm:ss or m:ss): 0:13.53
        Average shared text size (kbytes): 0
        Average unshared data size (kbytes): 0
        Average stack size (kbytes): 0
        Average total size (kbytes): 0
        Maximum resident set size (kbytes): 131952
        Average resident set size (kbytes): 0
        Major (requiring I/O) page faults: 0
        Minor (reclaiming a frame) page faults: 58974
        Voluntary context switches: 3
        Involuntary context switches: 26
        Swaps: 0
        File system inputs: 0
        File system outputs: 1968
        Socket messages sent: 0
        Socket messages received: 0
        Signals delivered: 0
        Page size (bytes): 4096
        Exit status: 0

Probably the most useful indicator above is Major (requiring I/O) page faults as this indicates whether the operating system is having to load pages of data from the disk because it no longer resides in RAM. This will cause a speed penalty.

In our example the code and data requirement is small so no page faults occur. If you have a memory-bound process, or several programs that use variable and large amounts of RAM, you might find that this gives you a clue as to which program is slowed down by disk accesses at the operating system level because parts of it have been swapped out of RAM to disk.

cProfile is a built-in profiling tool in the standard library. It hooks into the virtual machine in CPython to measure the time taken running every function that it sees. This introduces a greater overhead and you get more information as a result. Sometimes the additional information can lead to surprising insights into your code.

cProfile is one of three profilers in the standard library (the others are hotshot and profile). hotshot is experimental code, profile is the original pure-Python profiler. cProfile has the same interface as profile, is supported and is the default profiling tool. If you’re curious about the history of these libraries see Armin Rigo’s request to include cProfile in the standard library from 2005 ^[6].

A good practice when profiling is to generate a hypothesis about the speed of parts of your code before you profile it. Ian likes to print out the code snippet in question and to annotate it. Forming a hypothesis ahead of time means you can measure how wrong you are (you will get it wrong!) and improve your intuition about certain coding styles.

Always be driven by results that you have measured and always start with some quick-and-dirty profiling to make sure you’re addressing the right area. There’s nothing more humbling than cleverly optimizing a section of code only to realise (hours or days later) that you’d missed the slowest part of the process and haven’t really addressed the underlying problem at all.

Our hypothesis: We know that calculate_z_serial_purepython is likely to be the slowest part of the code. In that function we do a lot of dereferencing and make many calls to basic arithmetic operators and the abs function. These will probably show up as consumers of CPU resources.

Below we’ll use the cProfile module to run a variant of the code. The output is spartan but helps us figure out where to analyse further.

The -s cumulative flag tells cProfile to sort by cumulative time spent inside each function - this gives us a view into the slowest parts of a section of code. The cProfile output is written to screen directly after our usual print results.

$ python -m cProfile -s cumulative julia1_nopil.py
...
         36221992 function calls in 19.664 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.034    0.034   19.664   19.664 julia1_nopil.py:1(<module>)
        1    0.843    0.843   19.630   19.630 julia1_nopil.py:23(calc_pure_python)
        1   14.121   14.121   18.627   18.627 julia1_nopil.py:9
                                              (calculate_z_serial_purepython)
 34219980    4.487    0.000    4.487    0.000 {abs}
  2002000    0.150    0.000    0.150    0.000 {method 'append' of 'list' objects}
        1    0.019    0.019    0.019    0.019 {range}
        1    0.010    0.010    0.010    0.010 {sum}
        2    0.000    0.000    0.000    0.000 {time.time}
        4    0.000    0.000    0.000    0.000 {len}
        1    0.000    0.000    0.000    0.000 {method 'disable' of
                                               '_lsprof.Profiler' objects}

Sorting by cumulative time gives us an idea about where the majority of execution time is spent. The above result shows us that 36221992 function calls occurred in just over 19 seconds (this time includes the overhead of using cProfile). Previously our code took around 13 seconds to execute - we’ve just added a 5 second penalty by measuring how long each function takes to execute.

We can see that the entry point to the code julia1_cprofile.py on line 1 takes a total of 19 seconds, this is just the __main__ call to calc_pure_python. ncalls is 1, this line is only executed once.

Inside calc_pure_python the call to calculate_z_serial_purepython consumes 18.6 seconds, both functions are only called once. We can derive that approximately 1 second is spent on lines of code inside calc_pure_python separate to calling the CPU intensive calculate_z_serial_purepython function. We can’t derive which lines however take the time inside the function using cProfile.

Inside calculate_z_serial_purepython the time spent on lines of code (without calling other functions) is 14.1 seconds. This function makes 34,219,980 calls to abs which take a total of 4.4 seconds along with some other calls which do not cost much time.

What about the {abs} call? This line is measuring the individual calls to the abs function inside calculate_z_serial_purepython. Whilst the per-call cost is negligible (it is recorded as 0.000 seconds) the total time for 34,219,980 calls is 4.4 seconds. We couldn’t predict in advance exactly how many calls would be made to abs as the Julia function has non-predictable dynamics (that’s why it is so interesting to look at).

At best we could have said that it will be called a minimum of 1,000,000 times as we’re calculating 1000*1000 pixels. At most it would be 300,000,000 as we calculate 1,000,000 pixels with a maximum of 300 iterations. 34 million calls is roughly 10% of the worst case.

If we looked at the original greyscale image and in our mind’s eye we squashed the white parts together and into a corner, we can estimate that the expensive white region accounts for roughly 10% of the rest of the image.

The next line in the profiled output {method 'append' of 'list' objects} details the creation of 2,002,000 list items.

This creation of the 2,002,000 items is occurring in calc_pure_python during the setup phase.

The zs and cs will be 1000*1000 items each and these are built from a list of 1,000 x and 1,000 y co-ordinates, in total this is 2,002,000 calls to append.

It is important to note that this cProfile output is not ordered by parent functions, it is summarising the expense of all functions in the executed block of code. Figuring out what is happening on a line-by-line basis is very hard with cProfile as we only get profile information for the function calls themselves, not each line within the function.

Inside calculate_z_serial_purepython we can now account for {abs} and {range} and in total these two functions cost approximately 4.5 seconds. We know that calculate_z_serial_purepython costs 18.6 seconds.

The final line of the profiling output refers to lsprof - this is the original name of the tool that evolved into cProfile and can be ignored.

To get some more control over the results of cProfile we can write a statistics file and then analyse it in Python:

$ python -m cProfile -o profile.stats julia1.py

We can load this into Python with the following, it’ll give us the same cumulative time report as before:

In [1]: import pstats
In [2]: p = pstats.Stats("profile.stats")
In [3]: p.sort_stats("cumulative")
Out[3]: <pstats.Stats instance at 0x177dcf8>
In [4]: p.print_stats()
Tue Jan  7 21:00:56 2014    profile.stats

         36221992 function calls in 19.983 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.033    0.033   19.983   19.983 julia1_nopil.py:1(<module>)
        1    0.846    0.846   19.950   19.950 julia1_nopil.py:23(calc_pure_python)
        1   13.585   13.585   18.944   18.944 julia1_nopil.py:9(calculate_z_serial_purepython)
 34219980    5.340    0.000    5.340    0.000 {abs}
  2002000    0.150    0.000    0.150    0.000 {method 'append' of 'list' objects}
        1    0.019    0.019    0.019    0.019 {range}
        1    0.010    0.010    0.010    0.010 {sum}
        2    0.000    0.000    0.000    0.000 {time.time}
        4    0.000    0.000    0.000    0.000 {len}
        1    0.000    0.000    0.000    0.000 {method 'disable' of
                                               '_lsprof.Profiler' objects}

To trace which functions we’re profiling we can print the caller information. In the following two listings we can see that calculate_z_serial_purepython is the most expensive function and it is called from one place. If it was called from many places, these listings might help us narrow down the location of the most expensive parents:

In [5]: p.print_callers()
   Ordered by: cumulative time

Function                                was called by...
                                                 ncalls  tottime  cumtime
julia1_nopil.py:1(<module>)                   <-
julia1_nopil.py:23(calc_pure_python)          <-       1    0.846   19.950
                                                  julia1_nopil.py:1(<module>)
julia1_nopil.py:9(calculate_z_serial_purepython)  <-       1   13.585   18.944
                                                  julia1_nopil.py:23(calc_pure_python)
{abs}                                         <- 34219980    5.340    5.340
                                                  julia1_nopil.py:9(calculate_..._purepython)
{method 'append' of 'list' objects}           <- 2002000    0.150    0.150
                                                  julia1_nopil.py:23(calc_pure_python)
{range}                                       <-       1    0.019    0.019
                                                  julia1_nopil.py:9(calculate_..._purepython)
{sum}                                         <-       1    0.010    0.010
                                                  julia1_nopil.py:23(calc_pure_python)
{time.time}                                   <-       2    0.000    0.000
                                                  julia1_nopil.py:23(calc_pure_python)
{len}                                         <-       2    0.000    0.000
                                                  julia1_nopil.py:9(calculate_..._purepython)
                                                       2    0.000    0.000
                                                  julia1_nopil.py:23(calc_pure_python)
{method 'disable' of '_lsprof.Profiler' objects}  <-

We can flip this around the other way to show which functions call other functions:

In [6]: p.print_callees()
   Ordered by: cumulative time

Function                                          called...
                                                      ncalls  tottime  cumtime
julia1_nopil.py:1(<module>)                       ->       1    0.846   19.950
                                                      julia1_nopil.py:23(calc_pure_python)
julia1_nopil.py:23(calc_pure_python)              ->       1   13.585   18.944
                                                      julia1_nopil.py:9(calculate_..._purepython)
                                                           2    0.000    0.000
                                                      {len}
                                                     2002000    0.150    0.150
                                                      {method 'append' of 'list' objects}
                                                           1    0.010    0.010
                                                      {sum}
                                                           2    0.000    0.000
                                                      {time.time}
julia1_nopil.py:9(calculate_z_serial_purepython)  -> 34219980    5.340    5.340
                                                      {abs}
                                                           2    0.000    0.000
                                                      {len}
                                                           1    0.019    0.019
                                                      {range}
{abs}                                             ->
{method 'append' of 'list' objects}               ->
{range}                                           ->
{sum}                                             ->
{time.time}                                       ->
{len}                                             ->
{method 'disable' of '_lsprof.Profiler' objects}  ->

cProfile is rather verbose and you need a side screen to see it without lots of word-wrap. Since it is built-in it is a convenient tool for quickly identifying bottlenecks. Tools like line_profiler, heapy and memory_profiler which we discuss later in this chapter will then help you to drill-down to the specific lines that you should pay attention to.

runsnake is a visualisation tool for the profile statistics created by cProfile - you can quickly get a sense of which functions are most expensive just by looking at the diagram that’s generated.

Use runsnake to get a high level understanding of a cProfile statistics file especially when you’re investigating a new and large code-base, it’ll give you a feel for the areas that you should focus on. It might also reveal areas that you hadn’t realized would be expensive so it might highlight some quick wins for you to focus on.

You can also use it when discussing the slow areas of code in a team as it is easy to discuss the results.

Installation: pip install runsnake

Note that it requires wxPython and this can be a pain to install into a virtualenv. Ian has resorted to installing this globally on more than one occasion just to analyse a profile file, rather than to try to get it running in a virtualenv.

In Figure 2-5 we have the visual plot of the previous cProfile data. The visual inspection should make it easier to quickly understand that calculate_z_serial_purepython takes the majority of the time and that only a part of the execution time is due to calling other functions (the only one that is significant is abs). Visually you can see that there’s little point investing time in the setup routine as the vast majority of the execution time is in the calculation routine.

Figure 2-5. RunSnakeRun visualising a cProfile profile file.

With runsnake you can click into functions and drill into complex nested calls. When discussing the reasons for slow execution in a piece of code in a team this tool is invaluable.

In Ian’s opinion Robert Kern’s line_profiler is the strongest tool for identifying the cause of CPU-bound problems in Python code. It works by profiling individual functions on a line-by-line basis so you should start with cProfile and use the high-level view to guide which functions to profile with line_profiler.

It is worthwhile printing and annotating versions of the output from this tool as you modify your code so you have a record of changes (successful or not) that you can quickly refer to. Don’t rely on your memory when you’re working on line-by-line changes.

Installation: pip install line_profiler

A decorator (@profile) is used to mark the chosen function. The kernprof.py script is used to execute your code, the CPU time and other statistics for each line of the chosen function are recorded.

The arguments are -l for line-by-line (rather than function level) profiling and -v for a verbose output. Without -v you receive an .lprof output which you can later analyse with the line_profiler module. Below in Example 2-6 we’ll do a full run on our CPU-bound function.

$ kernprof.py -l -v julia1_lineprofiler.py
...
Wrote profile results to julia1_lineprofiler.py.lprof
Timer unit: 1e-06 s

File: julia1_lineprofiler.py
Function: calculate_z_serial_purepython at line 9
Total time: 100.81 s

Line #      Hits   Per Hit   % Time  Line Contents
==================================================
     9                               @profile
    10                               def calculate_z_serial_purepython(maxiter, zs, cs):
    11                                   """Calculate output list using
                                         Julia update rule"""
    12         1    6870.0      0.0      output = [0] * len(zs)
    13   1000001       0.8      0.8      for i in range(len(zs)):
    14   1000000       0.8      0.8          n = 0
    15   1000000       0.8      0.8          z = zs[i]
    16   1000000       0.8      0.8          c = cs[i]
    17  34219980       1.1     36.2          while abs(z) < 2 and n < maxiter:
    18  33219980       1.0     32.6              z = z * z + c
    19  33219980       0.8     27.2              n += 1
    20   1000000       0.9      0.9          output[i] = n
    21         1       4.0      0.0      return output

Introducing kernprof.py adds a substational time cost to the run time. In the above example calculate_z_serial_purepython takes 100 seconds, this is up from 13 seconds using simple print statements and 19 seconds using cProfile. The gain is that we get a line-by-line breakdown of where the time is spent inside the function.

The % Time column is the most helpful - we can see that 36% of the time is spent on the while testing. We don’t know whether the first statement (abs(z) < 2) is more expensive than the second (n < maxiter). Inside the loop we see that the update to z is also fairly expensive. Even n += 1 is expensive! Python’s dynamic lookup machinery is at work for every loop even though we’re using the same types for each variable in each loop - this is where compiling and type specialization (Chapter 7) give us a massive win. The creation of the output list and the updates on line 20 are relatively cheap compared to the cost of the while loop.

The obvious way to further analyse the while statement is to break it up. Whilst there has been some discussion in the Python community around the idea of re-writing the pyc files with more detailed information for multi-part single line statements, we are unaware of any production tools that offer a more fine-grained analysis than line_profiler.

Below in Example 2-7 we break the while logic into several statements. This additional complexity will increase the run time of the function as we have more lines of code to execute and it might help us to understand the costs incurred in this part of the code.

$ kernprof.py -l -v julia1_lineprofiler2.py
...
Wrote profile results to julia1_lineprofiler2.py.lprof
Timer unit: 1e-06 s

File: julia1_lineprofiler2.py
Function: calculate_z_serial_purepython at line 9
Total time: 184.739 s

Line #      Hits    Per Hit   % Time  Line Contents
===================================================
     9                                @profile
    10                                def calculate_z_serial_purepython(maxiter, zs, cs):
    11                                    """Calculate output list using
                                          Julia update rule"""
    12         1     6831.0      0.0      output = [0] * len(zs)
    13   1000001        0.8      0.4      for i in range(len(zs)):
    14   1000000        0.8      0.4          n = 0
    15   1000000        0.9      0.5          z = zs[i]
    16   1000000        0.8      0.4          c = cs[i]
    17  34219980        0.8     14.9          while True:
    18  34219980        1.0     19.0              not_yet_escaped = abs(z) < 2
    19  34219980        0.8     15.5              iterations_left = n < maxiter
    20  34219980        0.8     15.1              if not_yet_escaped
                                                      and iterations_left:
    21  33219980        1.0     17.5                  z = z * z + c
    22  33219980        0.9     15.3                  n += 1
    23                                            else:
    24   1000000        0.8      0.4                  break
    25   1000000        0.9      0.5          output[i] = n
    26         1        5.0      0.0      return output

This version takes 184 seconds to execute whilst the previous version took 100 seconds. Other factors did complicate the analysis. In this case having extra statements that have to be executed 34,219,980 times each slows down the code. If we hadn’t used kernprof.py to investigate the line-by-line effect of this change we might have drawn other conclusions about the reason for the slowdown as we’d have lacked the necessary evidence.

At this point it makes sense to step back to the earlier timeit technique to test the cost of individual expressions.

>>> z = 0+0j  # a point in the middle of our image
>>> %timeit abs(z) < 2  # tested inside IPython

10000000 loops, best of 3: 119 ns per loop

>>> n = 1
>>> maxiter = 300
>>>%timeit n < maxiter

10000000 loops, best of 3: 77 ns per loop

From this simple analysis it looks as though the logic test on n is almost two times faster than the call to abs. Since the order of evaluation for Python statements is both left-to-right and opportunistic it makes sense to put the cheapest test on the left side of the equation. On every 1 in 301 tests for each co-ordinate the n < maxiter test will be False so Python wouldn’t need to evaluate the other side of the and operator.

We never know whether abs(z) < 2 will be false until we evaluate it and our earlier observations for this region of the complex plane suggest it is True around 10% of the time for all 300 iterations. If we wanted to have a strong understanding of the time complexity of this part of the code it would make sense to continue the numerical analysis. In this situation however we want an easy check to see if we can get a quick win.

We can form a new hypothesis stating “by swapping the order of the operators in the while statement we will achieve a reliable speed-up”. We can test this hypothesis using kernprof.py but the additional overheads of profiling this way might add too much noise. Instead we can use an earlier version of the code, running a test comparing these two versions with while abs(z) < 2 and n < maxiter: against while n < maxiter and abs(z) < 2:.

The result was a fairly stable improvement of approximately 0.4 seconds. This result is obviously minor and is very problem specific, using a more suitable approach to solve this problem (e.g. swapping to use Cython or PyPy - see Chapter 7) would yield greater gains.

We can be confident in our result because:

For completeness we can run a final kernprof.py on the two main functions including our optimization to confirm that we have a complete picture of the overall complexity of our code:

$ kernprof.py -l -v julia1_lineprofiler3.py
...
Wrote profile results to julia1_lineprofiler3.py.lprof
Timer unit: 1e-06 s

File: julia1_lineprofiler3.py
Function: calculate_z_serial_purepython at line 9
Total time: 99.7097 s

Line #      Hits   Per Hit   % Time  Line Contents
==================================================
     9                               @profile
    10                               def calculate_z_serial_purepython(maxiter,
                                                                       zs, cs):
    11                                   """Calculate output list using
                                         Julia update rule"""
    12         1    6831.0      0.0      output = [0] * len(zs)
    13   1000001       0.8      0.8      for i in range(len(zs)):
    14   1000000       0.8      0.8          n = 0
    15   1000000       0.9      0.9          z = zs[i]
    16   1000000       0.8      0.8          c = cs[i]
    17  34219980       1.0     35.9          while n < maxiter and abs(z) < 2:
    18  33219980       1.0     32.0              z = z * z + c
    19  33219980       0.8     27.9              n += 1
    20   1000000       0.9      0.9          output[i] = n
    21         1       5.0      0.0      return output

Having swapped the two components of the while test on line 17 in Example 2-8 we see a modest improvement from 36.1% of the execution time to 35.9% (this result was stable over repeated runs).

File: julia1_lineprofiler3.py
Function: calc_pure_python at line 24
Total time: 195.218 s

Line #      Hits   Per Hit   % Time  Line Contents
==================================================
    24                               @profile
    25                               def calc_pure_python(draw_output,
                                                          desired_width, max_iterations):
...
    44         1       1.0      0.0      zs = []
    45         1       1.0      0.0      cs = []
    46      1001       1.1      0.0      for ycoord in y:
    47   1001000       1.1      0.5          for xcoord in x:
    48   1000000       1.5      0.8              zs.append(
                                                        complex(xcoord, ycoord))
    49   1000000       1.6      0.8              cs.append(
                                                        complex(c_real, c_imag))
    50
    51         1      51.0      0.0      print "Length of x:", len(x)
    52         1      11.0      0.0      print "Total elements:", len(zs)
    53         1       6.0      0.0      start_time = time.time()
    54         1  191031307.0     97.9   output =
                                           calculate_z_serial_purepython
                                           (max_iterations, zs, cs)
    55         1       4.0      0.0      end_time = time.time()
    56         1       2.0      0.0      secs = end_time - start_time
    57         1      58.0      0.0      print calculate_z_serial_purepython
                                           .func_name + " took", secs, "seconds"
    58
                                         # this sum is expected for 1000^2 grid with 300 iterations
    59         1    9799.0      0.0      assert sum(output) == 33219980

As expected in Example 2-9 calculate_z_serial_purepython takes most (97%) of the time of its parent function. The list creation steps are minor in comparison.

Just as Robert Kern’s line_profiler package measures CPU usage, the memory_profiler module by Fabian Pedregosa and Philippe Gervais measures memory usage on a line-by-line basis. By understanding the memory usage characteristics of your code you could ask yourself two questions:

memory_profiler operates in a very similar way to line_profiler but runs far more slowly. If you install the psutil package (optional but recommended) then memory_profiler will run faster. Memory profiling may easily make your code run 10 to 100 times slower. In practice you will probably use memory_profiler occasionally and line_profiler (for CPU profiling) more frequently.

Installation: pip install memory_profiler (and optionally pip install psutil)

As mentioned above the implementation of memory_profiler is not as performant as the implementation of line_profiler, it may make sense to run your tests on a smaller problem which completes in a useful amount of time. Overnight runs might be sensible for validation but you need quick and reasonable iterations to diagnose problems and hypothesise solutions.

The code in Example 2-10 uses the full 1000x1000 grid, the statistics took about 1.5 hours to collect on Ian’s laptop. When dealing with memory allocation you must be aware that the situation is not as clear-cut as it is with CPU usage.

Generally it is more efficient to over-allocate memory to a process into a local pool which can be used at leisure as memory allocation operations are relatively expensive. Garbage collection is not instant and so objects may be unavailable but still in the garbage collection pool for some time.

The outcome of this is that it is hard to really understand what is happening with memory usage and release inside a Python program as a line of code may not allocate a deterministic amount of memory as observed from outside the process. Observing the gross trend over a set of lines is likely to lead to better insight than observing the behaviour of just one line.

Inside calculate_z_serial_purepython on line 12 below we see that the allocation of 1,000,000 items causes approximately 7MB of RAM to be added to this process. This does not mean that the output list is definitely 7MB in size, just that the process grew by approximately 7MB during the internal allocation of the list. On line 13 we see that the process grew by approximately a further 32MB inside the loop, this can be attributed to the call to range. RAM tracking is further discussed in Example 11-1, the difference between 7MB and 32MB is due to the contents of the two lists.

In the parent process on line 46 we see that the allocation of the zs and cs lists consumes approximately 79MB. Again it is worth noting that this is not necessarily the true size of the arrays, just the size that the process grew by once these lists had been created.

$ python -m memory_profiler julia1_memoryprofiler.py
...
Line #    Mem usage    Increment   Line Contents
================================================
     9   89.934 MiB    0.000 MiB   @profile
    10                             def calculate_z_serial_purepython(maxiter, zs, cs):
    11                                 """Calculate output list using...
    12   97.566 MiB    7.633 MiB       output = [0] * len(zs)
    13  130.215 MiB   32.648 MiB       for i in range(len(zs)):
    14  130.215 MiB    0.000 MiB           n = 0
    15  130.215 MiB    0.000 MiB           z = zs[i]
    16  130.215 MiB    0.000 MiB           c = cs[i]
    17  130.215 MiB    0.000 MiB           while n < maxiter and abs(z) < 2:
    18  130.215 MiB    0.000 MiB               z = z * z + c
    19  130.215 MiB    0.000 MiB               n += 1
    20  130.215 MiB    0.000 MiB           output[i] = n
    21  122.582 MiB   -7.633 MiB       return output

Line #    Mem usage    Increment   Line Contents
================================================
    24   10.574 MiB -112.008 MiB   @profile
    25                             def calc_pure_python(draw_output,
                                                        desired_width,
                                                        max_iterations):
    26                                 """Create a list of complex ...
    27   10.574 MiB    0.000 MiB       x_step = (float(x2 - x1) / ...
    28   10.574 MiB    0.000 MiB       y_step = (float(y1 - y2) / ...
    29   10.574 MiB    0.000 MiB       x = []
    30   10.574 MiB    0.000 MiB       y = []
    31   10.574 MiB    0.000 MiB       ycoord = y2
    32   10.574 MiB    0.000 MiB       while ycoord > y1:
    33   10.574 MiB    0.000 MiB           y.append(ycoord)
    34   10.574 MiB    0.000 MiB           ycoord += y_step
    35   10.574 MiB    0.000 MiB       xcoord = x1
    36   10.582 MiB    0.008 MiB       while xcoord < x2:
    37   10.582 MiB    0.000 MiB           x.append(xcoord)
    38   10.582 MiB    0.000 MiB           xcoord += x_step
    ...
    44   10.582 MiB    0.000 MiB       zs = []
    45   10.582 MiB    0.000 MiB       cs = []
    46   89.926 MiB   79.344 MiB       for ycoord in y:
    47   89.926 MiB    0.000 MiB           for xcoord in x:
    48   89.926 MiB    0.000 MiB               zs.append(complex(xcoord, ycoord))
    49   89.926 MiB    0.000 MiB               cs.append(complex(c_real, c_imag))
    50
    51   89.934 MiB    0.008 MiB       print "Length of x:", len(x)
    52   89.934 MiB    0.000 MiB       print "Total elements:", len(zs)
    53   89.934 MiB    0.000 MiB       start_time = time.time()
    54                                 output = calculate_z_serial...
    55  122.582 MiB   32.648 MiB       end_time = time.time()
    ...

Another way to visualise the change in memory use is to sample over time and plot the result. memory_profiler has a utility called mprof, it is used once to sample the memory usage and a second time to visualise the samples. It samples by time and not by line so it barely impacts the run-time of the code.

Figure 2-6 is created using mprof run julia1_memoryprofiler.py, this writes a statistics file which is then visualised using mprof plot.

Our two functions are bracketed, this shows where in time they are entered and we can see the growth in RAM as they run. Inside calculate_z_serial_purepython we can see the steady increase in RAM usage throughout the execution of the function, this is caused by all the small objects (int and float types) that are created.

Figure 2-6. memory_profiler report using mprof.

In addition to observing the behaviour at the function level we can add labels using a context manager. The following snippet in Example 2-11 is used to generate the graph in Figure 2-7. The time.sleep(1) is an artificial addition to make the graph easier to understand. We can see the create_output_list label, it occurs momentarily after calculate_z_serial_purepython and the process has been allocated more RAM. We then pause for a second to make the graph clearer.

After the label create_range_of_zs we see a large and quick increase in RAM usage, in the modified code below you can see this label when we create the iterations list. Rather than using xrange instead we’ve used range - the diagram should make it clear that a large list of 1,000,000 elements is instantiated just for the purposes of generating an index and that this is an inefficient approach that will not scale to larger list sizes (we will run out of RAM!). The allocation of the memory used to hold this list will itself take a small amount of time which contributes nothing useful to this function.

@profile
def calculate_z_serial_purepython(maxiter, zs, cs):
    """Calculate output list using Julia update rule"""
    with profile.timestamp("create_output_list"):
        output = [0] * len(zs)
    time.sleep(1)
    with profile.timestamp("create_range_of_zs"):
        iterations = range(len(zs))
        with profile.timestamp("calculate_output"):
            for i in iterations:
                n = 0
                z = zs[i]
                c = cs[i]
                while n < maxiter and abs(z) < 2:
                    z = z * z + c
                    n += 1
                output[i] = n
    return output

In the calculate_output block which runs for most of the graph we see a very slow, linear increase in RAM usage, this will be from all of the temporary numbers used in the inner loops. Using the labels really helps us to understand at a fine-grained level where memory is being consumed.

Figure 2-7. memory_profiler report using mprof with labels.

Finally we could change the range call into an xrange and in Figure 2-8 we see the corresponding decrease in RAM usage in our inner loop.

Figure 2-8. memory_profiler report using mprof with labels.

If we want to measure the RAM used by several statements we can use the IPython magic %memit, it works just like %timeit. In Chapter 11 we look at using %memit to measure the memory cost of lists and discuss various ways of using RAM more efficiently.

The Guppy project has a heap inspection tool, it lets us look at the number and size of each object on Python’s heap. Looking inside the interpreter and understanding what’s held in memory is extremely useful in rare but difficult debugging sessions where you really need to know how many objects are in use and whether they get garbage collected at appropriate times. If you have a difficult memory leak (probably because references to your objects remain hidden in a complex system) then this is the tool that’ll get you to the root of the problem.

If you’re reviewing your code to see if it is generating as many objects as you predict then this tool is very useful - the results might surprise you and that could lead to new optimization opportunities.

Installation: pip install guppy

The listing in Example 2-12 is a slightly modified version of the Julia code, the heap object hpy is included in calc_pure_python and we print out the state of the heap at three places.

def calc_pure_python(draw_output, desired_width, max_iterations):
    ...
    while xcoord < x2:
        x.append(xcoord)
        xcoord += x_step

    from guppy import hpy; hp = hpy()
    print "heapy after creating y and x lists of floats"
    h = hp.heap()
    print h
    print

    zs = []
    cs = []
    for ycoord in y:
        for xcoord in x:
            zs.append(complex(xcoord, ycoord))
            cs.append(complex(c_real, c_imag))

    print "heapy after creating zs and cs using complex numbers"
    h = hp.heap()
    print h
    print

    print "Length of x:", len(x)
    print "Total elements:", len(zs)
    start_time = time.time()
    output = calculate_z_serial_purepython(max_iterations, zs, cs)
    end_time = time.time()
    secs = end_time - start_time
    print calculate_z_serial_purepython.func_name + " took", secs, "seconds"

    print
    print "heapy after calling calculate_z_serial_purepython"
    h = hp.heap()
    print h
    print

Example 2-13 shows that after creating the zs and cs lists the memory use is more interesting as it has grown by approximately 80MB due to 2,000,000 complex objects consuming 64,000,000 bytes. These complex numbers represent the majority of the memory usage at this stage. If you wanted to optimize the memory usage in this program then this result would be revealing - we can see both how many objects are being stored and their overall size.

guppy $ python julia1_guppy.py
heapy after creating y and x lists of floats
Partition of a set of 27293 objects. Total size = 3416032 bytes.
 Index  Count   %     Size   % Cumulative  % Kind (class / dict of class)
     0  10960  40  1050376  31   1050376  31 str
     1   5768  21   465016  14   1515392  44 tuple
     2    199   1   210856   6   1726248  51 dict of type
     3     72   0   206784   6   1933032  57 dict of module
     4   1592   6   203776   6   2136808  63 types.CodeType
     5    313   1   201304   6   2338112  68 dict (no owner)
     6   1557   6   186840   5   2524952  74 function
     7    199   1   177008   5   2701960  79 type
     8    124   0   135328   4   2837288  83 dict of class
     9   1045   4    83600   2   2920888  86 __builtin__.wrapper_descriptor
<91 more rows. Type e.g. '_.more' to view.>

heapy after creating zs and cs using complex numbers
Partition of a set of 2027301 objects. Total size = 83671256 bytes.
 Index  Count   %     Size   % Cumulative  % Kind (class / dict of class)
     0 2000000  99 64000000  76  64000000  76 complex
     1    185   0 16295368  19  80295368  96 list
     2  10962   1  1050504   1  81345872  97 str
     3   5767   0   464952   1  81810824  98 tuple
     4    199   0   210856   0  82021680  98 dict of type
     5     72   0   206784   0  82228464  98 dict of module
     6   1592   0   203776   0  82432240  99 types.CodeType
     7    319   0   202984   0  82635224  99 dict (no owner)
     8   1556   0   186720   0  82821944  99 function
     9    199   0   177008   0  82998952  99 type
<92 more rows. Type e.g. '_.more' to view.>

Length of x: 1000
Total elements: 1000000
calculate_z_serial_purepython took 13.2436609268 seconds

heapy after calling calculate_z_serial_purepython
Partition of a set of 2127696 objects. Total size = 94207376 bytes.
 Index  Count   %     Size   % Cumulative  % Kind (class / dict of class)
     0 2000000  94 64000000  68  64000000  68 complex
     1    186   0 24421904  26  88421904  94 list
     2 100965   5  2423160   3  90845064  96 int
     3  10962   1  1050504   1  91895568  98 str
     4   5767   0   464952   0  92360520  98 tuple
     5    199   0   210856   0  92571376  98 dict of type
     6     72   0   206784   0  92778160  98 dict of module
     7   1592   0   203776   0  92981936  99 types.CodeType
     8    319   0   202984   0  93184920  99 dict (no owner)
     9   1556   0   186720   0  93371640  99 function
<92 more rows. Type e.g. '_.more' to view.>

In the third section after calculating the Julia result we have used 94MB. In addition to the complex numbers we now have a large collection of integers in addition to more items stored in lists.

hpy.setrelheap() could be used to create a check-point of the memory configuration so subsequent calls to hpy.heap() will generate a delta from this check-point. This way you can avoid seeing the internals of Python and prior memory setup before the point of the program that you’re analysing.

Rober Brewer’s dowser hooks into the namespace of the running code and provides a real-time view of instantiated variables via a CherryPy interface in a web-browser. Each object that is being tracked has an associated sparkline graphic so you can watch to see if the quantities of certain objects are growing. This is useful for debugging long-running processes.

If you have a long-running process and you expect different memory behavior to occur depending on the actions you take in the application (e.g. with a web server you might upload data or cause complex queries to run) you can confirm this interactively, there’s an example in Figure 2-9.

Figure 2-9. Several sparklines shown through CherryPy with dowser.

To use it we add a convenience function, shown in Example 2-14, to the Julia code which can start the CherryPy server.

def launch_memory_usage_server(port=8080):
    import cherrypy
    import dowser

    cherrypy.tree.mount(dowser.Root())
    cherrypy.config.update({
        'environment': 'embedded',
        'server.socket_port': port
    })

    cherrypy.engine.start()

Before we begin our intensive calculations we launch the CherryPy server, shown in Example 2-15. Once we have completed our calculations we can keep the console open using the sleep, this leaves the CherryPy process running so we can continue to introspect the state of namespace.

 ...
        for xcoord in x:
            zs.append(complex(xcoord, ycoord))
            cs.append(complex(c_real, c_imag))

    launch_memory_usage_server()

...
    output = calculate_z_serial_purepython(max_iterations, zs, cs)
...
    print "now waiting..."
    while True:
        time.sleep(1)

By following the TRACE links in Figure 2-10 we can view the contents of each list object. We could further drill down into each list, this is like using an interactive debugger in an IDE but you can do this on a deployed server without an IDE.

Figure 2-10. 1,000,000 items in a list with dowser.

So far we’ve reviewed various ways to measure the cost of Python code (both for CPU and RAM usage). We haven’t yet looked at the underlying bytecode used by the virtual machine, understanding what’s happening “under the hood” helps to build a mental model of what’s happening in slow functions and it’ll help when you come to compile your code. Let’s introduce some bytecode.

The dis module lets us inspect the underlying bytecode that we run inside the stack based CPython virtual machine. Having an understanding of what’s happening in the virtual machine that runs your higher level Python code will help you to understand why some styles of coding are faster than others. It will also help when you come to use a tool like Cython which steps outside of Python and generates C code.

The dis module is built-in. You can pass it code or a module and it will print out a disassembly. In Example 2-16 we disassemble the outer loop of our CPU-bound function.

In [1]: import dis
In [2]: import julia1_nopil

In [3]: dis.dis(julia1_nopil.calculate_z_serial_purepython)
 11           0 LOAD_CONST               1 (0)
              3 BUILD_LIST               1
              6 LOAD_GLOBAL              0 (len)
              9 LOAD_FAST                1 (zs)
             12 CALL_FUNCTION            1
             15 BINARY_MULTIPLY
             16 STORE_FAST               3 (output)

 12          19 SETUP_LOOP             123 (to 145)
             22 LOAD_GLOBAL              1 (range)
             25 LOAD_GLOBAL              0 (len)
             28 LOAD_FAST                1 (zs)
             31 CALL_FUNCTION            1
             34 CALL_FUNCTION            1
             37 GET_ITER
        >>   38 FOR_ITER               103 (to 144)
             41 STORE_FAST               4 (i)

 13          44 LOAD_CONST               1 (0)
             47 STORE_FAST               5 (n)

# ...
# We'll snip the rest of the inner loop for brevity!
# ...

 19     >>  131 LOAD_FAST                5 (n)
            134 LOAD_FAST                3 (output)
            137 LOAD_FAST                4 (i)
            140 STORE_SUBSCR
            141 JUMP_ABSOLUTE           38
        >>  144 POP_BLOCK

 20     >>  145 LOAD_FAST                3 (output)
            148 RETURN_VALUE

The output of the above code is fairly straightforward, if terse. The first column contains line numbers that relate to our original file. The second column contains several >> symbols, these are the destinations for jump points elsewhere in the code. The third column is the operation address along with the operation name. The fourth column contains the parameters for the operation. The fifth column contains annotations to help line up the bytecodes with the original Python parameters.

Refer back to Example 2-3 to match the bytecode to the corresponding Python code. The bytecode starts by putting the constant value 0 onto the stack and then it builds a single element list. Next it searches the namespaces to find the len function, puts it on the stack, searches the namespaces again to find zs and then puts that onto the stack. On line 12 it calls the len function from the stack which consumes the zs reference in the stack, then it applies a binary multiply to the last two arguments (the length of zs and the single element list) and stores the result in output. That’s the first line of our Python function now dealt with. Follow the next block of bytecode to understand the behavior of the second line of Python code (the outer for loop).

Having introduced bytecode we can now ask - what’s the bytecode and time cost of writing an function out explicitly versus using built-ins to perform the same task?

def fn_expressive(upper = 1000000):
    total = 0
    for n in xrange(upper):
        total += n
    return total

def fn_terse(upper = 1000000):
    return sum(xrange(upper))

print "Functions return the same result:", fn_expressive() == fn_terse()

Functions return the same result:
True

%timeit fn_expressive()
%timeit fn_terse()

10 loops, best of 3: 42 ms per loop
100 loops, best of 3: 12.3 ms per loop

import dis
print fn_expressive.func_name
dis.dis(fn_expressive)

fn_expressive
  2           0 LOAD_CONST               1 (0)
              3 STORE_FAST               1 (total)

  3           6 SETUP_LOOP              30 (to 39)
              9 LOAD_GLOBAL              0 (xrange)
             12 LOAD_FAST                0 (upper)
             15 CALL_FUNCTION            1
             18 GET_ITER
        >>   19 FOR_ITER                16 (to 38)
             22 STORE_FAST               2 (n)

  4          25 LOAD_FAST                1 (total)
             28 LOAD_FAST                2 (n)
             31 INPLACE_ADD
             32 STORE_FAST               1 (total)
             35 JUMP_ABSOLUTE           19
        >>   38 POP_BLOCK

  5     >>   39 LOAD_FAST                1 (total)
             42 RETURN_VALUE

print fn_terse.func_name
dis.dis(fn_terse)

fn_terse
  8           0 LOAD_GLOBAL              0 (sum)
              3 LOAD_GLOBAL              1 (xrange)
              6 LOAD_FAST                0 (upper)
              9 CALL_FUNCTION            1
             12 CALL_FUNCTION            1
             15 RETURN_VALUE

If you aren’t already unit-testing your code then you are probably hurting your longer-term productivity. Ian (blushing) is embarrassed to note that once he spent a day optimising his code, having disabled unit-tests because they were inconvenient, only to discover that his significant speed-up result was due to breaking a part of the algorithm he was improving. You do not need to make this mistake even once.

In addition to unit-testing you should also strongly consider coverage.py. It checks to see which lines of code are exercised by your tests and identifies the sections that have no coverage. This quickly lets you figure out whether you’re testing the code that you’re about to optimize, such that any mistakes that might creep in during the optimisation process are quickly caught.

# ex.py
import unittest

@profile
def some_fn(nbr):
    return nbr * 2


class TestCase(unittest.TestCase):
    def test(self):
        result = some_fn(2)
        self.assertEquals(result, 4)

$ nosetests ex.py
E
======================================================================
ERROR: Failure: NameError (name 'profile' is not defined)
...
NameError: name 'profile' is not defined

Ran 1 test in 0.001s

FAILED (errors=1)

# for line_profiler
if '__builtin__' not in dir() or not hasattr(__builtin__, 'profile'):
    def profile(func):
        def inner(*args, **kwargs):
            return func(*args, **kwargs)
        return inner

$ kernprof.py -v -l ex.py
Line #      Hits         Time  Per %%HTMLit   % Time  Line Contents
==============================================================
    11                                           @profile
    12                                           def some_fn(nbr):
    13         1            3      3.0    100.0      return nbr * 2


$ nosetests ex.py
.
Ran 1 test in 0.000s

# for memory_profiler
if 'profile' not in dir():
    def profile(func):
        def inner(*args, **kwargs):
            return func(*args, **kwargs)
        return inner

 $ python -m memory_profiler ex.py
...
Line #    Mem usage    Increment   Line Contents
================================================
    11   10.809 MiB    0.000 MiB   @profile
    12                             def some_fn(nbr):
    13   10.809 MiB    0.000 MiB       return nbr * 2


$ nosetests ex.py
.
Ran 1 test in 0.000s

Profiling requires some time and concentration. You will stand a better chance of understanding your code if you separate the section you want to test from the main body of your code. By separating it you can unit-test the code to preserve correctness and you can pass in realistic fabricated data to exercise the inefficiencies you want to address.

Do remember to disable any BIOS-based accelerators, they will only confuse your results. On Ian’s laptop the Intel TurboBoost feature can temporarily accelerate a CPU above its normal maximum speed if it is cool enough. This means that a cool CPU may run the same block of code faster than a hot CPU. Your Operating System may also control the clock speed - a laptop on battery power is likely to more aggressively control CPU speed than a laptop on mains power. To create a more stable benchmarking configuration we:

Try to hypothesise the expected behaviour of your code and then validate (or disprove!) the hypothesis with the result of a profiling step. Your choices will not change (you can only use the profiled results to drive your decisions) but your intuitive understanding of the code will improve and this will pay off in future projects as you are more likely to make performant decisions. Of course you will verify these performant decisions by profiling as you go.

Do not skimp on the preparation. If you try to performance test code deep inside a larger project without separating it from the larger project then you are likely to witness side effects which will side-track your efforts. It is likely to be harder to unit-test a larger project when you’re making fine-grained changes and this may further hamper your efforts. Side effects could include other threads and processes impacting CPU, memory, network and disk activity which will skew your results.

For web servers investigate dowser and dozer, you can use these to visualise in real-time the behavior of objects in the namespace. Definitely consider separating the code you want to test out of the main web application if possible, it will make profiling significantly easier.

Make sure your unit-tests exercise all the code paths in the code that you’re analysing. Anything you don’t test but is used in your benchmarking may cause subtle errors that will slow down your progress. Use coverage.py to confirm that your tests are covering all the code paths.

Unit-testing a complicated section of code that generates a large numerical output may be difficult. Do not be afraid to output a text file of results to run through diff or to use a pickled object. For numeric optimization problems Ian likes to create long text files of floating point numbers and to use diff - minor rounding errors show up immediately and can be rare in the output.

If your code might be subject to numerical rounding issues due to subtle changes then you are better off with a large output that can be used for a before-and-after comparison. One cause of rounding errors is the difference in floating point precision between CPU registers and main memory. Running your code through a different code path can cause subtle rounding errors which might later confound you - it is better to be aware of this as soon as they occur.

Obviously it makes sense to use a source code control tool whilst you a profiling and optimizing. Branching is cheap and it will preserve your sanity.

Having looked at profiling techniques you should have all the tools you need to identify bottlenecks around CPU and RAM usage in your code. Next we’ll look at how Python implements the most common containers so you can make sensible decisions about representing larger collections of data.

As alluded to before, we can achieve better search performance if we first sort our data so that all elements to the left of a particular item are smaller than it (or larger). The comparison is done through the __eq__ and __lt__ magic function of the object and can be user defined if using custom objects.

The two elements necessary are the sorting algorithm and then the searching algorithm. Python lists have a built in sorting algorithm that uses Tim Sort. Tim Sort can sort through a list in O(n) in the best case (and O(n logn) in the worst case). It achieves this performance by utilizing multiple types of sorting algorithms and using heuristics to guess which algorithm will perform the best given the data (more specifically, it hybridizes insertion and merge sort algorithms).

Once a list has been sorted, we can use a binary search algorithm (Example 3-3) to find our desired element, which has an average case complexity of O(logn). It achieves this by first looking at the middle of the list and comparing this value with the desired value. If this midpoint’s value is less than our desired value, then we consider the left half of the list and continue halving the list like this until the value is found. As a result, we do not need to read all values in the list as was necessary for the linear search and instead only read a small subset of them.

def binary_search(needle, haystack):
    imin, imax = 0, len(haystack)
    while True:
        if imin >= imax:
            return -1
        midpoint = (imin + imax) // 2
        if haystack[midpoint] > needle:
            imax = midpoint
        elif haystack[midpoint] < needle:
            imin = midpoint+1
        else:
            return midpoint

This method allows us to find elements in a list without resorting to the potentially heavy weight solution of a dictionary. Especially when the list of data that is being operated on is intrinsically sorted, it is more efficient to simply do a binary search on the list to find an object (and get the O(logn) complexity of the search) rather than converting your data to a dictionary and then doing a lookup on it (although a dictionary lookup takes O(1), converting to a dictionary takes O(n) and a dictionary’s restriction of no repeating keys may be undesirable).

In addition, the bisect module simplifies much of this process by giving easy methods to add elements into a list while maintaining its sorting, in addition to finding elements using a heavily optimized binary search. This is done by providing alternative functions to append that, instead of adding the new element to the end of the list, adds the element into the correct sorted placement. With the list always being sorted, we can easily find the elements we are looking for (examples of this can be found in the documentation for the bisect module ^[8]). In addition, we can use bisect to find the closest element to what we are looking for very quickly (Example 3-4). This can be extremely useful when comparing two datasets that are similar but do not have exactly the same results.

import bisect
import random

def find_closest(haystack, needle):
    # bisect.bisect_left will return the first value in the haystack that is greater
    # than or equal to the needle
    i = bisect.bisect_left(haystack, needle)
    if i == len(haystack):
        return i - 1
    elif haystack[i] == needle:
        return i
    elif i > 0:
        j = i - 1
        # since we know value is larger than needle (and vice versa for the
        # value at j), we don't need to use absolute values here
        if haystack[i] - needle > needle - haystack[j]:
            return j
    return i

important_numbers = []
for i in xrange(10):
    new_number = random.randint(0, 1000)
    bisect.insort(important_numbers, new_number)

# important_numbers will be already in order because we inserted new elements
# with bisect.insort
print important_numbers

closest_index = find_closest(important_numbers, -250)
print "Closest value to -250: ", important_numbers[closest_index]

closest_index = find_closest(important_numbers, 500)
print "Closest value to 500: ", important_numbers[closest_index]

closest_index = find_closest(important_numbers, 1100)
print "Closest value to 1100: ", important_numbers[closest_index]

In general this touches on a fundamental rule of writing efficient code: pick the right data structure and stick with it! Although there may be more efficient data structures for particular operations, the cost of the conversion to that data structure may defeat any efficiency boost.

If lists and tuples both use the same underlying data structure, what are the differences between the two? Summarized, the main differences are:

These differences instill the philosophical difference between the two: tuples are for describing multiple properties of one unchanging thing and lists can be used to store collections of data about completely disparate objects. For example, the parts of a telephone number are perfect for a tuple: they won’t change and if they do they represent a new object or a different phone number. Similarly, the coefficients of a polynomial fit a tuple since different coefficients represent a different polynomial. On the other hand, the names of the people currently reading this book are better suited for a list: the data is constantly changing both in content and in size but is still always representing the same idea.

It is important to note that both lists and tuples can take mixed types. This can, as we will see, introduce quite some overhead and reduce some potential optimizations. These overheads can be removed if we also force all the data in our list to be of the same type. In Chapter 6 we will talk about reducing both the memory and computational overhead by using numpy. In addition, other packages such as blist and array can also reduce these overheads for other, non-numerical, situations. This alludes to a major point in performant programming that we will touch on in later chapters: abstract code will be much slower than code specifically designed to solve a particular problem.

In addition, the immutability of a tuple, opposed to a list which can be resized and changed, makes them a very lightweight data structure. This means that there isn’t much overhead in memory when storing them and operations with them are quite straightforward. With lists, as we will learn, their mutability [Link to Come] comes at a price of extra memory needed to store them and extra computations necessary when using them.

Lists as dynamic arrays

Once we create a list, we are free to change its contents as needed,

>>> numbers = [5, 8, 1, 3, 2, 6]
>>> numbers[2] = 2*numbers[0]  # 
>>> numbers
[5, 8, 10, 3, 2, 6]

As described before, this operation is O(1) because we can find the data stored within the 0th and 2nd elements immediately.

In addition, we can append new data to a list and grow it’s size:

>>> len(numbers)
6
>>> numbers.append(42)
>>> numbers
[5, 8, 10, 3, 2, 6, 42]
>>> len(numbers)
7

This is possible because dynamic arrays support a resize operation which increases the capacity of the array. When a list is of size N is first appended to, python must create a new list that is big enough to hold the original N items in addition to the extra one that is being appended. However, instead of allocating N+1 items, M items are actually allocated where M>N in order to provide extra headroom for future appends. Then, the data from the old list is copied to the new list and the old list is destroyed. The philosophy is that one append is probably the beginning of many appends and by requesting extra space we can reduce the number of times this allocation must happen and the total number of memory copies that are necessary. This is quite important since memory copies can be quite expensive, especially when list sizes start growing. Figure Figure 3-2 shows what this overallocation looks like in python2.7. The formula dictating this growth is:

Example 3-5. List overallocation equation in py2.7

M = (N >> 3) + (N < 9 ? 3 : 6)

N	0	1-4	5-8	9-16	17-25	26-35	36-46	…	991-1120
M	0	4	8	16	25	35	46	…	1120

Figure 3-2. This graph shows how many extra elements are being allocated to a list of a particular size

As we append data, we utilize the extra space and increase the effective size of the list, N. As a result, N grows as we append new data until N == M. At this point, there is no extra space to insert new data into and we must create a new list with more extra space. This new list has extra headroom as given by the equation above and we copy the old data into the new space.

This sequence of events is shown visually in Figure Figure 3-3 below. The figure follows the various operations being performed on list l in Example 3-6.

Example 3-6. Resizing of a list

l = [1, 2]
for i in range(3, 7):
    l.append(i)

Figure 3-3. Example of how a list is mutated on multiple appends

Note

This extra allocation happens on the first append. When a list is directly created, as in the example above, only the number of elements needed are allocated.

While the amount of extra headroom allocated is generally quite small, it can add up. This effect becomes especially pronounced when maintaining many small lists or when keeping a particularly large list. If we are storing 1,000,000 lists, each containing 10 elements, we would suppose that 10,000,000 elements worth of memory is being used. In actuality, however, up to 16,000,000 elements could have been allocated if the append operator was used to construct the list. Similarly, for a large list of 100,000,000 elements, we actually have 112,500,007 elements allocated!

>>> t = (1,2,3,4)
>>> t[0] = 5
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'tuple' object does not support item assignment

>>> t1 = (1,2,3,4)
>>> t2 = (5,6,7,8)
>>> t1 + t2
(1, 2, 3, 4, 5, 6, 7, 8)

>>> %timeit l = [0,1,2,3,4,5,6,7,8,9]
1000000 loops, best of 3: 285 ns per loop
>>> %timeit t = (0,1,2,3,4,5,6,7,8,9)
10000000 loops, best of 3: 55.7 ns per loop

Lists and tuples are fast and low overhead objects to use when your data already has an intrinsic ordering to it. This intrinsic ordering allows you to sidestep the search problem in these structures, if the ordering is known beforehand then lookups are O(1) instead of needing an expensive O(N) linear search. While lists can be resized, care needs to be taken to properly understand how much overallocation is happening to ensure that the dataset can still fit in memory. On the other hand, tuples can be created quickly and without the added overhead of lists, at the cost of not being modifiable.

In the next chapter, we go over the computational properties of dictionaries that solve the search/lookup problems with un-ordered data at the cost of overhead.

Dictionaries and sets use a hash table in order to achieve its O(1) lookup and insertions. This efficiency is a result of a very clever usage of hash functions [Link to Come] in order to create list-like indexes for data. By turning the data’s key into something that can be used like a list index, we can get the same the same performance as a list. In addition, instead of having to refer to data by a numerical index, which itself implies some ordering to the data, we can refer to it by this arbitrary key.

Inserting and Retrieving

In order to create a hash table from scratch, we start with some allocated memory, similar to what we started with for arrays. For an array, if we wanted to insert data, we simply find the smallest unused bucket and insert our data there (and resize if necessary). For hash tables, we must first figure out the placement of the data in this contiguous chunk of memory.

The placement of the new data is contingent on two properties of the data you are inserting: the hashed value of the key and how the value compares to other objects. This is because when we insert data, first the key is hashed and masked so that it turns into an effective index in an array.^[10] The mask makes sure that the hash value, which can take the value of any integer, fits within the allocated number of buckets. So, if we have allocated 8 blocks of memory and our hash value is 28975, we consider the bucket at index 28975 & 0b111 = 7. If, however, our dictionary has grown to have 256 items then the mask becomes 0b11111111. Now we must check if this bucket is already in use. If it is empty, we can insert the key and the value into this block of memory. We store the key so that we can make sure we are retrieving the correct value on lookups. If it is in use and the value of the bucket is equal to the value we wish to insert (a comparison done with the cmp builtin), then the key/value pair is already in the hash table and we can return. However, if the values don’t match then we must find a new place to put the data.

In order to find a new index, we compute a new index using a simple linear function, a method called probing. In addition, python’s probing mechanism adds a contribution from higher order bits of the original hash (recall that for a table of length 8 we only considered the last 3 bits of the hash for the initial index through the use of a mask value of mask = 0b111 = bin(8 - 1)). Using these higher order bits gives each hash a different sequence of next possible hashes, which helps to avoid future collisions. There is a lot of freedom when picking the algorithm to generate a new index, however it is quite important that the scheme visits every possible index in order to evenly distribute the data in the table. How well distributed the data is throughout the hash table is called the “load factor” and is related to the entropy of the hash function. The following pseudo-code illustrates the calculation of hash indexes used in cpython2.7:

Example 4-4. Dictionary lookup sequence

def index_sequence(key, mask=0b111, PERTURB_SHIFT=5):
    perturb = hash(key) 
    i = perturb & mask
    yield i
    while True:
        i = ((i << 2) + i + perturb + 1)
        perturb >>= PERTURB_SHIFT
        yield i & mask

hash returns an integer while the actual C code in CPython uses an unsigned integer. Because of this, the above pseudocode doesn’t 100% replicate the behaviour in CPython, however it is a good approximation.

This probing is a modification of the naive method of “linear probing”. In linear probing, we simply yield the values i = (5*i + 1) & mask where i is initialized to the hash value of the key and the value 5 is unimportant ^[11]. An important thing to note is that linear probing only deals with the several bytes of the hash and disregards the rest (ie: for a dictionary with 8 elements, we only look at the last 3 bits). This means that if two item’s hash have the same last 3 binary digits, we will not only have a collision, but the sequence of probed indexes will be the same. The perturbed scheme that python uses will start taking into consideration more bits from the items hash in order to resolve this problem.

A similar procedure is done when we are performing lookups on a specific key: the given key is transformed into an index and that index is looked at. If the key in that index matches (recall that we also store the original key when doing insert operations), then we can return that value. If it doesn’t, we keep creating new indexes using the same scheme as above until we either find the data or hit an empty bucket. If we hit an empty bucket, we can conclude that the data does not exist in the table.

Figure 4-1, illustrates the process of adding some data into a hash table. Here, we chose the hashing function to be simply the ascii value of the first letter of the key using Python’s ord function which returns an integer ordinal for a single character string. As we’ll see later in Hash functions and Entropy, python provides hashing functions for most of its types so it is normally not necessary to provide one yourself except in extreme situations. Insertion of the key “Barcelona” causes a collision and a new index is calculated using the scheme the scheme in Example 4-4. This dictionary can also be created in python using the following code:

Example 4-5. Custom hashing function

class City(str):
    def __hash__(self):
        return ord(self[0])

# We create a dictionary where we assign arbitrary values to cities.
data =  {
    City("Rome"): 4,
    City("San Francisco"): 3,
    City("New York"): 5,
    City("Barcelona"): 2,
}

In this case “Barcelona” and “Rome” cause the hash collision (Figure 4-1 shows the outcome of this insertion). We see this because for a dictionary with 4 elements we have a mask value of 0b111. As a result, “Barcelona” will try to use index ord("B") & 0b111 = 66 & 0b111 = 0b1000010 & 0b111 = 0b010 = 2. Similarly “Rome” will try to use the index ord("R") & 0b111 = 82 & 0b111 = 0b1010010 & 0b111 = 0b010 = 2.

Dictionary Mechanics

Finding an element. Using the dictionary created in Example 4-5, what would a lookup on the key “Johannesburg” look like? What indices would be checked?

Deleting an element. Using the dictionary created in Example 4-5, how would you handle the deletion of the key “Rome”? How would subsequent lookups for the key “Rome” or “Barcelona” be handled?

Hash Collisions. Considering the dictionary created in Example 4-5, how many hash collisions could you expect if 500 cities, all starting with an upper-case letter, were added into a hash table? How about 1000 cities? Can you think of a way of lowering this number?

Discussion. For 500 cities there would be approximately 474 dictionary elements which collided with a previous value (500 - 26) with each hash having 500 / 26 = 19.2 cities associated with it. For 1000 cities, 974 elements would collide and each hash would have 1000 / 26 = 38.4 cities associated with it. This is because the hash is simply based on the numerical value of the first letter which can only take a value from A-Z allowing for only 26 independent hash values. This means that a lookup in this table could require as many as 38 subsequent lookups to find the correct value. In order to fix this, we must increase the number of possible hash values by considering other aspects of the city in the hash. The default hash function on a string considers every character in order to maximize the number of possible values. See Hash functions and Entropy for more explanation.

Figure 4-1. The resulting hash table from inserting with collisions

16, 64, 256, 1024, 4096, 16384, 65536, 131072, 262144, ...

class Point(object):
    def __init__(self, x, y):
        self.x, self.y = x, y

>>> p1 = Point(1,1)
>>> p2 = Point(1,1)
>>> set([p1, p2])
set([<__main__.Point at 0x1099bfc90>, <__main__.Point at 0x1099bfbd0>])
>>> Point(1,1) in set([p1, p2])
False

class Point(object):
    def __init__(self, x, y):
        self.x, self.y = x, y

    def __hash__(self):
        return hash((self.x, self.y))

    def __eq__(self, other):
        return self.x == other.x and self.y == other.y

>>> p1 = Point(1,1)
>>> p2 = Point(1,1)
>>> set([p1, p2])
set([<__main__.Point at 0x109b95910>])
>>> Point(1,1) in set([p1, p2])
True

def twoletter_hash(key):
    offset = ord('a')
    k1, k2 = key
    return (ord(k2) - offset) + 26 * (ord(k1) - offset)

import string
import timeit

class BadHash(str):
    def __hash__(self):
        return 42

class GoodHash(str):
    def __hash__(self):
        """
        This is a slightly optimized version of twoletter_hash
        """
        return ord(self[1]) + 26 * ord(self[0]) - 2619

baddict = set()
gooddict = set()
for i in string.ascii_lowercase:
    for j in string.ascii_lowercase:
        key = i + j
        baddict.add(BadHash(key))
        gooddict.add(GoodHash(key))

badtime = timeit.repeat(
    "key in baddict",
    setup = "from __main__ import baddict, BadHash; key = BadHash('zz')",
    repeat = 3,
    number = 1000000,
)
goodtime = timeit.repeat(
    "key in gooddict",
    setup = "from __main__ import gooddict, GoodHash; key = GoodHash('zz')",
    repeat = 3,
    number = 1000000,
)

print "Min lookup time for baddict: ", min(badtime)
print "Min lookup time for gooddict: ", min(goodtime)

# Results:
#   Min lookup time for baddict:  16.3375990391
#   Min lookup time for gooddict:  0.748275995255

Doing a lookup on a dictionary is fast however unnecessarily doing it will slow down your code, just as any extraneous lines will. One way this surfaces itself is in Python’s namespace management which heavily uses dictionaries to do its lookups.

Whenever a variable, function or module is being invoked in Python, there is a hierarchy where these objects are looked for. First, Python looks inside of the locals() array which has entries for all local variables. Python works hard to make local variable lookups fast and this is the only part of the chain that doesn’t require a dictionary lookup. If it doesn’t exist there, then the globals() dictionary is searched. Finally, if we don’t find it there, the __builtin__ object is searched. It is important to note that while locals() and globals() are explicitly dictionaries, while __builtin__ is technically a module object, when searching __builtin__ for a given property we are just doing a dictionary lookup inside of its locals() map (this is the case for all module object and class objects!).

To make this clearer, let’s look at a simple example of calling functions that are defined in different scopes (Example 4-8). We can disassemble the functions with the dis module (Example 4-9) to get a better understanding of how these namespace lookups are happening.

import math
from math import sin

def test1(x):
    """
    >>> %timeit test1(123456)
    1000000 loops, best of 3: 381 ns per loop
    """
    return math.sin(x)

def test2(x):
    """
    >>> %timeit test2(123456)
    1000000 loops, best of 3: 311 ns per loop
    """
    return sin(x)

def test3(x, sin=math.sin):
    """
    >>> %timeit test3(123456)
    1000000 loops, best of 3: 306 ns per loop
    """
    return sin(x)

>>> dis.dis(test1)
  9           0 LOAD_GLOBAL              0 (math)  // Dictionary Lookup
              3 LOAD_ATTR                1 (sin)   // Dictionary Lookup
              6 LOAD_FAST                0 (x)     // Local Lookup
              9 CALL_FUNCTION            1
             12 RETURN_VALUE

>>> dis.dis(test2)
 15           0 LOAD_GLOBAL              0 (sin)  // Dictionary Lookup
              3 LOAD_FAST                0 (x)    // Local Lookup
              6 CALL_FUNCTION            1
              9 RETURN_VALUE

>>> dis.dis(test3)
 21           0 LOAD_FAST                1 (sin)  // Local Lookup
              3 LOAD_FAST                0 (x)    // Local Lookup
              6 CALL_FUNCTION            1
              9 RETURN_VALUE

The first function, test1, makes the call to sin by explicitly looking at the math library. This is also evident in the byte code that is produced, first a reference to the math module must be loaded, then we do an attribute lookup on this module until we finally have a reference to the sin function. This is done through 2 dictionary lookups, one to find the math module and one to find the sin function within the module.

On the other hand, test2 explicitly imports the sin function from the math module, then the function is directly accessible within the global namespace. This means we can avoid the lookup of the math module and the subsequent attribute lookup. However, we still must find the sin function within the global namespace. This is yet another reason to be explicit about what functions you are importing from a module. This practice not only makes code more readable because the reader knows exactly what functionality is required from external sources, but it also speeds up code!

Finally, test3 defines the sin function as a keyword argument with its default value being a reference to the sin function within the math module. While we still do need to find a reference to this function within the module, it only must happen when the test3 function is first defined. After this, the reference to the sin function is stored within the function definition as a local variable in the form of a default keyword argument. As mentioned before, local variables do not need a dictionary lookup to be found, they are stored in a very slim array which has very fast lookup times. Because of this, finding the function is quite fast!

While these effects are an interesting result of the way namespaces in python are managed, test3 is definitely not “Pythonic”. Luckily, these extra dictionary lookups only start to degrade performance when they are called a lot (in the innermost block of a very fast loop such as in the Julia set example). With this in mind, a more readable solution would be to set a local variable with the global reference before the loop is started. This still means we must do the global lookup once whenever the function is called however all the calls to that function in the loop will be made faster. This speaks to the fact that even minute slowdowns in code can be amplified if that code is being run millions of times. Even though a dictionary lookup may only take a several hundred nanoseconds, if we are looping millions of times over this lookup it can quickly add up. In fact, looking at Example 4-10 we see a 9.4% speedup simply by making the sin function local a tight loop which calls it.

from math import sin

def tight_loop_slow(iterations):
    """
    >>> %timeit tight_loop_slow(10000000)
    1 loops, best of 3: 2.21 s per loop
    """
    result = 0
    for i in xrange(iterations):
        # this call to sin requires a global lookup
        result += sin(i)

def tight_loop_fast(iterations):
    """
    >>> %timeit tight_loop_fast(10000000)
    1 loops, best of 3: 2.02 s per loop
    """
    result = 0
    local_sin = sin
    for i in xrange(iterations):
        # this call to local_sin requires a local lookup
        result += local_sin(i)

Dictionaries and sets provide a fantastic way to store data that can be indexed by a key. The way this key is used, through the hashing function, can greatly effect the resulting performance of the data structure. Furthermore, understanding how dictionaries work not only gives you a better understanding of how to organize your data, but also how to organize your code since dictionaries are a very intrinsic part of python’s internal functionality.

In the next chapter, we will explore generators which allows us to provide data to code in any dynamic order we want and without necessarily storing it in memory beforehand. This let’s us side-step many of the possible hurdles when using any of python’s intrinsic data structures.

Since we only need to store some version of a state and emit only the current value, generators are ideal for infinite series. The Fibonacci series is a great example—it is an infinite series with 2 state variables (the two last Fibonacci numbers).

def fibonacci():
    i, j = 0, 1
    while True:
        yield j
        i, j = j, i + j

We see here that although j is the value being emitted, we keep track of i as well since this holds the state of the Fibonacci series. The amount of state necessary for a calculation is quite important for generators because it translates into the actual memory footprint of the object. This makes it so that if we had a function that uses a lot of state and outputs a very small amount of data, it may be better to have this function pre-compute the list of data rather than have a generator for it.

One reason why generators aren’t used as much as they could be, is that a lot of the logic within them can be encapsulated in your logic code. This means that generators are really a way of organizing your code and having smarter loops. For example, we could answer the question “how many Fibonacci numbers below 5000 are odd” in multiple ways,

def fibonacci_naive():
    i, j = 0, 1
    count = 0
    while j <= 5000:
        if j % 2:
            count += 1
        i, j = j, i + j
    return count

def fibonacci_generator():
    count = 0
    for f in fibonacci():
        if f > 5000:
            break
        if f % 2:
            count += 1
    return count

from itertools import islice
def fibonacci_succinct():
    is_odd = lambda x : x % 2
    first_5000 = islice(fibonacci(), 0, 5000)
    return sum(1 for x in first_5000 if is_odd(x))

All of these methods have similar runtime properties (namely they all have the same memory footprint and the same performance), however the fibonacci_genarator function benefits from several things. Firstly, it is much more verbose than fibonacci_succinct which means it will be easy for another developer to debug and understand. This function mainly stands as a warning for the next section where we cover some common work flows using itertools—while the module greatly simplifies many simple actions with iterators, it can also quickly make Python code very un-Pythonic. Conversely, fibonacci_naive is doing multiple things at a time which hides the actual calculation it is doing! While it is obvious in the generator function we are iterating over the Fibonacci numbers, we are not over encumbered by the actual calculation. Lastly, fibonacci_generator is more generalizable. This function could be renamed num_odd_under_5000 and take in the generator by argument and thus work over any series.

One last benefit of the fibonacci_generator function is that it supports the notion that in computation there are two phases: generating data and transforming data. This function is very clearly a transformation on data while the fibonacci function generates it. This clear demarcation adds extra clarity and extra functionality by being able to move a transformative function to work on a new set of data, or adding new transformations onto existing data. This paradigm has always been important when creating complex programs, however generators facilitate this clearly by making generators responsible for creating the data, and normal functions responsible for acting on the generated data.

As touched on before, the way we get the memory benefits with a generator is by forcing us only to deal with the current value of interest. At any point in our calculation with a generator, we only have the current value and cannot reference any other items in the sequence (algorithms that perform this way are generally called “single pass” or “online”). This can make generators sometimes more difficult to use, however there are many modules and functions that can help.

The main library of interest is itertools in the standard library. It supplies generator versions of Python’s builtin functions map, reduce, filter and zip (called imap, ireduce, filter and izip in itertools) in addition to many other useful functions. Of particular note are:

Let’s build up an example of using generators in order to analyse a large dataset. Let’s say we had an analysis routine going over temporal data, one piece of data per second for the last 20 years—631,152,000 data points! The data is stored in a file, one second per line, and we cannot load the entire data-set into memory. If we wanted to do some simple anomaly detection, we could use generators and never allocate any lists!

The problem will be, given a datafile of the form “timestamp, value”, find days with a value 3 sigma away from that day’s mean. We start by writing the code that will read the file, line by line, and output that line’s value as python objects. We will also create a read_fake_data generator to generate fake data we can test our algorithms with. For this function we still take the argument filename so as to have the same function signature as read_data, however we will simply disregard it. These two functions are indeed lazily evaluated—we only read the next line in the file, or generate new fake data, when the next() property on the generator is called.

from random import normalvariate, rand
from itertools import count

def read_data(filename):
    with open(filename) as fd:
        for line in fd:
            data = line.strip().split(',')
            yield map(int, data)

def read_fake_data(filename):
    for i in count():
        sigma = rand() * 10
        yield (i, normalvariate(0, sigma))

Now, we can use the groupby function in itertools to group together timestamps that occur on the same day. This function works by taking in a sequence of items and a key used to group these items. The output is a generator which produces tuples whose items are the key for the group and a generator for the items in the group. As our key function, we will create a lambda function that returns a date object. These date objects are equal when they occur on the same day, which will group them by day! This “key” function could be anything—we could group our data by hour, by year, or by some property in the actual value. The only limitation is that groups will only be formed for data that is sequential. So, if we had the input be A A A A B B A A and had groupby group by the letter, we would get three groups, (A, [A, A, A, A]), (B, [B, B]), (A, [A, A]).

from datetime import date
from itertools import groupby

def day_grouper(iterable):
    key = lambda (timestamp, value) : date.fromtimestamp(timestamp)
    return groupby(iterable, key)

Now to do the actual anomaly detection. We will do this by iterating through a day’s values and keeping track of the mean and maximum values. The mean will be calculated using an online mean and standard deviation algorithm ^[13]. The maximum is kept because it is the best candidate for being anomalous—if the max is 3 sigma larger than the mean then we will will return the date object representing the day we just analysed. Otherwise, we will return False for posterity, however we could equally have just ended the function (and implicitly returned None). We output these values because this check_anomaly function is defined as a data filter—a function that returns a value that evaluates to True for data that should be kept and a value that evaluates to False if a value should be disregarded. This will allow us to filter through the original dataset and only keep days which pass our condition.

import math

def check_anomaly((day, day_data)):
    # We find the mean, standard deviation and maximum values for the day.
    # Using a single pass mean/standard deviation algorithm allows us to only
    # read through the day's data once.
    n = 0
    mean = 0
    M2 = 0
    max_value = None
    for timestamp, value in day_data:
        n += 1
        delta = value - mean
        mean = mean + delta/n
        M2 += delta*(value - mean)
        max_value = max(max_value, value)
    variance = M2/(n - 1)
    standard_deviation = math.sqrt(variance)

    # Here is the actual check of whether that day's data is anomalous.  If it
    # is, we return the value of the day, otherwise we return false
    if max_value > mean + 3 * standard_deviation:
        return day
    return False

One aspect of this function which may seem strange is the extra parenthesis in the parameter definition. This is not a typo but a result of this function taking input from the groupby generator. Recall that groupby returns a tuple, each one of which will become the parameters to this check_anomaly function. As a result, we must do tuple expansion in order to properly extract the group key and the group data. Since we are using ifilter, another way of dealing with this instead of having the tuple expansion inside the function definition is to define istarfilter which is similar to what istarmap does to imap (see the itertools documentation for more information).

Finally, we can put together the chain of generators to get the days that had anomalous data,

from itertools import ifilter, imap

data = read_data(data_filename)
data_day = day_grouper(data)
anomalous_dates = ifilter(None, imap(check_anomaly, data_day)) #

first_anomalous_date, first_anomalous_data = anomalous_dates.next()
print "The first anomalous date is: ", first_anomalous_date

This method very simply allows us to get the list of days that are anomalous without having to load the entire dataset. One thing to note is that the above code does not actually do any calculation, it simply sets up the pipeline to do the calculation. The file will never be read until we do anomalous_dates.next() or somehow iterate on the anomalous_dates generator. In fact, we only ever do the analysis when a new value is requested from anomalous_dates. So, if our full dataset has 5 anomalous dates, but in our code we retrieve one and then stop requesting new values, the file will only be read up to the point where this day’s data appears. This is called lazy evaluation—only the calculations that are explicitly requested are performed which can drastically reduce overall runtime if there is a early termination condition.

Another nicety about organizing analysis this way is it easily allows us to do more expansive calculations without having to rework large parts of code. For example, if we wanted to have a moving window of one day instead of chunking up by days, we can simply write a new day_grouper that does so and use it instead,

from datetime import datetime

def rolling_window_grouper(data, window_size=3600):
    window = tuple(islice(data, 0, window_size))
    while True:
        current_datetime = datetime.fromtimestamp(window[0][0])
        yield (current_datetime, window)
        window = window[1:] + (data.next(),)

Now, we simply replace the call to day_grouper in Example 5-5 with a call to rolling_window_grouper and we get the desired result. In this version we also see very explicitly the memory guarantee of this and the previous method—it will store only the window’s worth of data as state (in both cases one day or 3600 data points). This can be changed by opening the file multiple time and using this to point to the parts of the data we like (or using the linecache module), however this is only necessary if this sub-sampling of the dataset still doesn’t fit in memory.

A final note, in the rolling_window_grouper function, we perform many pop and append operations on the list window. This can be greatly optimized by using the deque object in the collections module. This object gives us O(1) removals and appends from the beginning or end of a list (while normal lists have O(1) for appends or removals from the end of the list and O(N) for the same operations on the beginning of the list). Using the deque object, we can append the new data to the right (or end) of the list and deque.popleft() to delete the leftmost (or beginning) of the list without having to allocate more space or perform long O(N) operations.

By formulating our anomaly finding algorithm as iterators, we are able to process much more data than could fit into memory. Not only are we able to process more data than can fit in RAM, but we also are able to do it faster than if we had used lists since we were able to avoid all the costly append operations.

Since iterators are a primitive type in Python, this should always be a go-to method when trying to reduce the memory footprint of an application. The benefits are that results are lazily evaluated, so you only ever process the data you need, and memory is saved since we don’t store previous results unless explicitly required to. In Chapter 11, we will talk about other methods that can be used for more specific problems and hopefully teach new ways of looking at problems when RAM is an issue.

Another benefit of thinking about solving problems using iterators is it prepares your code to be used on multiple CPU’s or multiple computers, as we will see in Chapter 9 and Chapter 10. As we discussed in Iterators for Infinite Series, when working with iterators you must always think about the various states that are necessary for your algorithm to work. Once you figure out how to package the state necessary for the algorithm to run, it doesn’t matter where it runs. We can see this sort of paradigm, for example, with the multiprocessing and ipython modules both using a map like function to launch parallel tasks.

In order to explore matrix and vector computation in this chapter, we will repeatedly use the example of diffusion in fluids. Diffusion is one of the mechanisms that moves fluids and tries to make it be uniformly mixed.

In this section we will explore the mathematics behind the diffusion equation. This may seem complicated, but don’t worry! We will quickly simplify this to make it more understandable. Also, it is important to note that while having a basic understanding of the final equation we are solving will be useful while reading this chapter, it is not completely necessary; the subsequent chapters will focus mainly on various formulations of the code, not the equation. However, having an understanding of the equations helps gain intuition about ways of optimizing our code. This is true in general—understanding the motivation behind your code and the intricacies of the algorithm will help you gain an intuition about methods of optimization.

One simple example of diffusion is dye in water: if you put several drops of dye into water at room temperature it will slowly move out until it fully mixes with the water. Since we are not stirring the water, nor is it warm enough to create convection currents, diffusion will be the main process mixing the two liquids. When solving these equations numerically, we pick what we want the initial conditions to look like and are able to evolve the initial condition forward in time to see what it will look like at a later time (see Figure 6-2).

All this being said, the most important thing to know about diffusion for our purposes is its formulation. Stated as a partial differential equation in 1 dimension it is written as,

Equation 6-1. Diffusion Equation

In this formulation, u is the vector representing the quantities we are diffusing. For example, we could have a vector with values of 0 where there is only water and 1 where there is only dye (and values in-between where there is mixing). In general, this will be a 2D or 3D matrix representing an actual area or volume of fluid. In this way, we could have u be a 3D matrix representing the fluid in a glass and instead of simply doing the second derivative along the x direction, we’d have to take it over all axes. In addition, D is a physical value that represents properties of the fluid we are simulating. A large value of D represents a fluid that can diffuse very easily. For simplicity, we will set D = 1 for our code, but still include it in the calculations.

What we will do is take the diffusion equation above, which is continuous in space and time, and approximate it to discrete volumes and discrete times. This is done using Euler’s method. Euler’s method simply take’s the derivative and writes it as a difference, such that

where dt is now a fixed number. This fixed number represents the time step size, or the resolution in time that we wish to solve this equation. It can be thought of as the frame rate of the movie we are trying to make. As the frame rate goes up (or dt goes down), we get a clearer picture of what happens. In fact, as dt approaches zero, Euler’s approximation becomes exact. We can thus rewrite the above equation to figure out what u(x, t+dt) is given u(x,t). What this means for us is that we can start with some initial state (u(x,0) representing the glass of water just as we add a drop of dye into it) and churn through the above mechanisms to “evolve” that initial state and see what it will look like at future times (u(x,dt)). This type of problem is called an “initial value problem” or “Cauchey Problem”.

Doing a similar trick for the derivative in x using the finite differences approximation, we arrive at the final equation,

Here, similar to how dt represents the frame rate, dx represents the resolution of the images—the smaller dx is the smaller a region every cell in our matrix represents. For simplicity, we will simply set D = 1 and dx = 1. These two values become very important when doing proper physical simulations, however since we are simply solving the diffusion equation for illustrative purposes they are not important to us.

Using this equation, we can solve almost any diffusion problem. However, there are some considerations we must make from this equation. Firstly, we said before that the spacial index in u (ie: the x parameter) will be represented as the indices into a matrix. What happens when we try to find the value at x-dx when x is at the beginning of the matrix? This problem is called the boundary condition. You can have fixed boundary conditions that say: any value out of the bounds of my matrix will be set to 0 (or any other value). Alternatively, you can have periodic boundary conditions that say that the values will wrap around (ie: if your matrix has length N, the value at index -1 is the same as at N-1, and the value at N+1 is the same at index 1. In order words, if you are trying to access the value at index i, you will get the value at index (i%N)).

Another consideration is how we are going to store the multiple time components of u. We could have one matrix for every time value we do our calculation for. At minimum, it seems that we will need two matrices; one matrix for the current state of the fluid and one for the next state of the fluid. As we’ll see, there are very drastic performance considerations for this particular question.

So, what does it look like to solve this problem in practice? Below is some pseudo code to illustrate the way we can use this equation to solve the problem

# Create the initial conditions
u = vector of length N
for i in range(N):
    u = 0 if there is water, 1 if there is dye

# Evolve the initial conditions
D = 1
t = 0
dt = 0.0001
while True:
    print "Current time is: %f" % t
    unew = vector of size N

    # update step for every cell
    for i in range(N):
        unew[i] = u[i] + dt * (u[(i+1)%N] + u[(i-1)%N] - 2 * u[i])
    # move the updated solution into u
    u = unew

    visualize(u)

This code will take some initial condition of the dye in water and tell us what the system looks like at every 0.0001 second interval in the future. The results of this can be seen in Figure 6-1 where we evolve a very concentrated drop of dye (represented by the top-hat function) into the future. We can see how far into the future, the dye becomes well mixed to the point where everywhere has a similar concentration of the dye.

Figure 6-1. Example of 1D diffusion

For the purposes of this chapter, we will be solving the 2D version of the above equation. All this means is that instead of operating over a vector (or in other words, a matrix with one index), we will be operating over a 2D matrix. The only change to the equations (and thus the subsequent code), is that we must now also take the second derivative in the y direction. This simply means that the original equation we are working with becomes,

Equation 6-2. Numerical diffusion equation in 2D

Which translates to the following pseudo code using the same methods we used before,

for i in range(N):
    for j in range(M):
        unew[i][j] = u[i][j] + dt * (                       \
            (u[(i+1)%N][j] + u[(i-1)%N][j] - 2 * u[i][j]) + \ # d^2 u / dx^2
            (u[i][(j+1)%M] + u[j][(j-1)%M] - 2 * u[i][j])   \ # d^2 u / dy^2
        )

We can put all of this together and write the full python 2D diffusion code which we will use as the basis for our benchmarks for the rest of this chapter. While the code looks more complicated, the results are similar to that of the 1D diffusion (as can be seen in Figure 6-2).

Figure 6-2. Example of diffusion for two sets of initial conditions

Let’s take python code from Example 6-1 with the 2D extension and break it up into chunks so that we can better analyse its runtime performance. The first step is to write out the evolution function that takes in the matrix and returns its evolved state,

grid_shape = (1024, 1024)

def evolve(grid, dt, D=1.0):
    xmax, ymax = grid_shape
    new_grid = [[0.0 for x in xrange(ymax)] for x in xrange(xmax)]
    for i in xrange(xmax):
        for j in xrange(ymax):
            grid_xx = grid[(i+1)%xmax][j] + grid[(i-1)%xmax][j] - 2.0 * grid[i][j]
            grid_yy = grid[i][(j+1)%ymax] + grid[i][(j-1)%ymax] - 2.0 * grid[i][j]
            new_grid[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt
    return grid

The global variable grid_shape designates how big of a region we be simulating and, as explained in Introduction to the Problem, we are using periodic boundary conditions (which is why we use modulo for the indices). In order to actually use this code, we must initialize a grid and call evolve on it. The following code is a very generic initialization procedure and will be reused throughout the chapter (its performance characteristics will not be analysed since it must only run once, as opposed to the evolve function, which is called repeatedly).

def run_experiment(num_iterations):
    # setting up initial conditions 
    grid = [[0.0 for x in xrange(grid_shape[1])] for x in xrange(grid_shape[0])]

    block_low = int(grid_shape[0] * .4)
    block_high = int(grid_shape[0] * .5)
    for i in xrange(block_low, block_high):
        for j in xrange(block_low, block_high):
            grid[i][j] = 0.005

    # Evolve the initial conditions
    start = time.time()
    for i in range(num_iterations):
        grid = evolve(grid, 0.1)
    return time.time() - start

The values for dt and grid elements have been chosen to be sufficiently small so that the algorithm is stable. See “Numerical Recipes” for a more in-depth treatment of this algorithm’s convergence characteristics.

$ kernprof.py -lv diffusion_python.py
Wrote profile results to diffusion_python.py.lprof
Timer unit: 1e-06 s

File: diffusion_python.py
Function: evolve at line 8
Total time: 16.1398 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     8                                           @profile
     9                                           def evolve(grid, dt, D=1.0):
    10        10           39      3.9      0.0      xmax, ymax = grid_shape
    11   2626570      2159628      0.8     13.4      new_grid = [[0.0 for x in xrange(ymax)] for x in xrange(xmax)]
    12      5130         4167      0.8      0.0      for i in xrange(xmax):
    13   2626560      2126592      0.8     13.2          for j in xrange(ymax):
    14   2621440      4259164      1.6     26.4              grid_xx = grid[(i+1)%xmax][j] + grid[(i-1)%xmax][j] - 2.0 * grid[i][j]
    15   2621440      4196964      1.6     26.0              grid_yy = grid[i][(j+1)%ymax] + grid[i][(j-1)%ymax] - 2.0 * grid[i][j]
    16   2621440      3393273      1.3     21.0              new_grid[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt
    17        10           10      1.0      0.0      return grid

from math import sin

def loop_slow(num_iterations):
    """
    >>> %timeit loop_slow(int(1e4))
    100 loops, best of 3: 2.67 ms per loop
    """
    result = 0
    for i in xrange(num_iterations):
        result += i * sin(num_iterations) #
    return result

def loop_fast(num_iterations):
    """
    >>> %timeit loop_fast(int(1e4))
    1000 loops, best of 3: 1.38 ms per loop
    """
    result = 0
    factor = sin(num_iterations)
    for i in xrange(num_iterations):
        result += i * factor
    return result

def evolve(grid, dt, out, D=1.0):
    xmax, ymax = grid_shape
    for i in xrange(xmax):
        for j in xrange(ymax):
            grid_xx = grid[(i+1)%xmax][j] + grid[(i-1)%xmax][j] - 2.0 * grid[i][j]
            grid_yy = grid[i][(j+1)%ymax] + grid[i][(j-1)%ymax] - 2.0 * grid[i][j]
            out[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt

def run_experiment(num_iterations):
    # setting up initial conditions
    scratch = [[0.0 for x in xrange(grid_shape[1])] for x in xrange(grid_shape[0])]
    grid = [[0.0 for x in xrange(grid_shape[1])] for x in xrange(grid_shape[0])]

    block_low = int(grid_shape[0] * .4)
    block_high = int(grid_shape[0] * .5)
    for i in xrange(block_low, block_high):
        for j in xrange(block_low, block_high):
            grid[i][j] = 0.005

    start = time.time()
    for i in range(num_iterations):
        evolve(grid, 0.1, scratch)
        grid, scratch = scratch, grid
    return time.time() - start

$ kernprof.py -lv diffusion_python_memory.py
Wrote profile results to diffusion_python_memory.py.lprof
Timer unit: 1e-06 s

File: diffusion_python_memory.py
Function: evolve at line 8
Total time: 13.3209 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     8                                           @profile
     9                                           def evolve(grid, dt, out, D=1.0):
    10        10           15      1.5      0.0      xmax, ymax = grid_shape
    11      5130         3853      0.8      0.0      for i in xrange(xmax):
    12   2626560      1942976      0.7     14.6          for j in xrange(ymax):
    13   2621440      4059998      1.5     30.5              grid_xx = grid[(i+1)%xmax][j]
                                                                       + grid[(i-1)%xmax][j]
                                                                       - 2.0 * grid[i][j]
    14   2621440      4038560      1.5     30.3              grid_yy = grid[i][(j+1)%ymax]
                                                                       + grid[i][(j-1)%ymax]
                                                                       - 2.0 * grid[i][j]
    15   2621440      3275454      1.2     24.6              out[i][j] = grid[i][j] + D
                                                                         * (grid_xx + grid_yy)
                                                                         * dt

The python code we wrote at Example 6-6 still has a problem that is at the heart of using python for these sorts of vectorized operations: python doesn’t natively support vectorization. This happens because of two things: python lists store pointers to the actual data and python bytecode is not optimized for vectorization so for loops cannot predict when using vectorization would be beneficial.

The fact that python lists store pointers means that instead of actually holding the data we care about, lists store a location where the data can actually be found. For most uses, this is good because it allows us to store whatever type of data we like inside of a list. However, when it comes to vector and matrix operations this is a source for a lot of performance degradation.

The degradation occurs because every time we want to fetch an element from the grid matrix, we must do multiple lookups. For example, doing grid[5][2] requires us to first do a list lookup for index 5 on the list grid. This will return a pointer to where the data at that location is stored. Then, we do another list lookup on this returned object for the element at index 2. Once we have this reference, we have the location where the actual data is stored.

The overhead for one such lookup is not big and can be, in most cases, disregarded. However, if the data we wanted was located in one contiguous block in memory we could move ALL of the data in one operation instead of needing 2 operations for each element. This is one of the major points with data fragmentation: when your data is fragmented you must move each piece over individually instead of moving the entire block over. This means you are invoking more memory transfer overhead and you are forcing the CPU to wait while data is being transfered. We will see with perf just how important this is when looking at the cache-misses.

This problem of getting the right data to the CPU when it is needed is related to the “Von Neumann bottleneck”. This refers to the fact that there is a limited bandwidth between the memory and the CPU as a result of the tiered memory architecture which modern computers use. If we could move data infinitely fast, we would not need any cache since the CPU could retrieve any data it needed instantly. This would be a state where the bottleneck is non-existent.

In order to deal with this, we pre-fetch data from the RAM and fill smaller but faster CPU caches with the data so that hopefully, when the CPU needs a piece of data, it will be located in a location that can be read from quickly. While this is a severely idealized way of looking at the architecture, we can still see some of the problems with it—how do we know what data will be needed in the future? The CPU does a good job with a mechanisms called branch prediction and pipelining which try to predict the next instruction and load the relevant memory into cache while still working on the current instruction. However, the best way to minimize the effects of the bottleneck is to be smart about how we allocate our memory and how we do our calculations over our data.

Probing how well memory is being moved to the CPU can be quite hard, however in Linux the perf tool can be used to get amazing amounts of insight regarding how the CPU is dealing with the program being run. For example, we can run perf on the pure python code from Example 6-6 and see just how efficiently the CPU is running our code,

$ perf stat -e cycles,stalled-cycles-frontend,stalled-cycles-backend,instructions,\
    cache-references,cache-misses,branches,branch-misses,task-clock,faults,minor-faults,\
    cs,migrations -r 3 python diffusion_python_memory.py

 Performance counter stats for 'python diffusion_python_memory.py' (3 runs):

   329,155,359,015 cycles                    #    3.477 GHz                      ( +-  0.47% )
    76,800,457,550 stalled-cycles-frontend   #   23.33% frontend cycles idle     ( +-  1.75% )
    46,556,100,820 stalled-cycles-backend    #   14.14% backend  cycles idle     ( +-  0.77% )
   598,135,111,009 instructions              #    1.82  insns per cycle
                                             #    0.13  stalled cycles per insn  ( +-  0.00% )
        35,497,196 cache-references          #    0.375 M/sec                    ( +-  0.94% )
        10,716,972 cache-misses              #   30.191 % of all cache refs      ( +- 15.51% )
   133,881,241,254 branches                  # 1414.067 M/sec                    ( +-  0.00% )
     2,891,093,384 branch-misses             #    2.16% of all branches          ( +-  0.12% )
      94678.127621 task-clock                #    0.999 CPUs utilized            ( +-  0.48% )
             5,439 page-faults               #    0.057 K/sec                    ( +-  0.01% )
             5,439 minor-faults              #    0.057 K/sec                    ( +-  0.01% )
               125 context-switches          #    0.001 K/sec                    ( +-  2.55% )
                 6 CPU-migrations            #    0.000 K/sec                    ( +- 16.67% )

      94.749389121 seconds time elapsed                                          ( +-  0.48% )

from array import array
import numpy

def norm_square_list(vector):
    """
    >>> vector = range(1000000)
    >>> %timeit norm_square_list(vector_list)
    1000 loops, best of 3: 1.16 ms per loop
    """
    norm = 0
    for v in vector:
        norm += v*v
    return norm

def norm_square_list_comprehension(vector):
    """
    >>> vector = range(1000000)
    >>> %timeit norm_square_list(vector_list)
    1000 loops, best of 3: 913 µs per loop
    """
    return sum([v*v for v in vector])

def norm_square_array(vector):
    """
    >>> vector = array('l', range(1000000))
    >>> %timeit norm_square_array(vector_array)
    1000 loops, best of 3: 1.44 ms per loop
    """
    norm = 0
    for v in vector:
        norm += v*v
    return norm

def norm_square_numpy(vector):
    """
    >>> vector = numpy.arange(1000000)
    >>> %timeit norm_square_numpy(vector_numpy)
    10000 loops, best of 3: 30.9 µs per loop
    """
    return numpy.sum(vector * vector) # 

def norm_square_numpy_dot(vector):
    """
    >>> vector = numpy.arange(1000000)
    >>> %timeit norm_square_numpy_dot(vector_numpy)
    10000 loops, best of 3: 21.8 µs per loop
    """
    return numpy.dot(vector, vector) #